{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# _*Pricing European Call Options*_ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "<br>\n",
    "Suppose a <a href=\"http://www.theoptionsguide.com/call-option.aspx\">European call option</a> with strike price $K$ and an underlying asset whose spot price at maturity $S_T$ follows a given random distribution.\n",
    "The corresponding payoff function is defined as:\n",
    "\n",
    "$$\\max\\{S_T - K, 0\\}$$\n",
    "\n",
    "In the following, a quantum algorithm based on amplitude estimation is used to estimate the expected payoff, i.e., the fair price before discounting, for the option:\n",
    "\n",
    "$$\\mathbb{E}\\left[ \\max\\{S_T - K, 0\\} \\right]$$\n",
    "\n",
    "as well as the corresponding $\\Delta$, i.e., the derivative of the option price with respect to the spot price, defined as:\n",
    "\n",
    "$$\n",
    "\\Delta = \\mathbb{P}\\left[S_T \\geq K\\right]\n",
    "$$\n",
    "\n",
    "The approximation of the objective function and a general introduction to option pricing and risk analysis on quantum computers are given in the following papers:\n",
    "\n",
    "- <a href=\"https://www.nature.com/articles/s41534-019-0130-6\">Quantum Risk Analysis. Woerner, Egger. 2018.</a>\n",
    "- <a href=\"https://quantum-journal.org/papers/q-2020-07-06-291/\">Option Pricing using Quantum Computers. Stamatopoulos et al. 2019.</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:09.974086Z",
     "start_time": "2020-07-13T23:35:09.967567Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "\n",
    "from qiskit import Aer, QuantumCircuit\n",
    "from qiskit.aqua.algorithms import IterativeAmplitudeEstimation\n",
    "from qiskit.circuit.library import LogNormalDistribution, LinearAmplitudeFunction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Uncertainty Model\n",
    "\n",
    "We construct a circuit factory to load a log-normal random distribution into a quantum state.\n",
    "The distribution is truncated to a given interval $[\\text{low}, \\text{high}]$ and discretized using $2^n$ grid points, where $n$ denotes the number of qubits used.\n",
    "The unitary operator corresponding to the circuit factory implements the following: \n",
    "\n",
    "$$\\big|0\\rangle_{n} \\mapsto \\big|\\psi\\rangle_{n} = \\sum_{i=0}^{2^n-1} \\sqrt{p_i}\\big|i\\rangle_{n},$$\n",
    "\n",
    "where $p_i$ denote the probabilities corresponding to the truncated and discretized distribution and where $i$ is mapped to the right interval using the affine map:\n",
    "\n",
    "$$ \\{0, \\ldots, 2^n-1\\} \\ni i \\mapsto \\frac{\\text{high} - \\text{low}}{2^n - 1} * i + \\text{low} \\in [\\text{low}, \\text{high}].$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:10.294975Z",
     "start_time": "2020-07-13T23:35:10.283925Z"
    }
   },
   "outputs": [],
   "source": [
    "# number of qubits to represent the uncertainty\n",
    "num_uncertainty_qubits = 3\n",
    "\n",
    "# parameters for considered random distribution\n",
    "S = 2.0       # initial spot price\n",
    "vol = 0.4     # volatility of 40%\n",
    "r = 0.05      # annual interest rate of 4%\n",
    "T = 40 / 365  # 40 days to maturity\n",
    "\n",
    "# resulting parameters for log-normal distribution\n",
    "mu = ((r - 0.5 * vol**2) * T + np.log(S))\n",
    "sigma = vol * np.sqrt(T)\n",
    "mean = np.exp(mu + sigma**2/2)\n",
    "variance = (np.exp(sigma**2) - 1) * np.exp(2*mu + sigma**2)\n",
    "stddev = np.sqrt(variance)\n",
    "\n",
    "# lowest and highest value considered for the spot price; in between, an equidistant discretization is considered.\n",
    "low  = np.maximum(0, mean - 3*stddev)\n",
    "high = mean + 3*stddev\n",
    "\n",
    "# construct A operator for QAE for the payoff function by\n",
    "# composing the uncertainty model and the objective\n",
    "uncertainty_model = LogNormalDistribution(num_uncertainty_qubits, mu=mu, sigma=sigma**2, bounds=(low, high))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEyCAYAAADOV2anAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3debgcVZnH8e9PEAgEMLKERSSACiJRMEGMoCSArDpBEKKiYxCNjAKOIgwiQgAdAWVxYBQjakRGoyKiIIhhSRDZ40JYomwhBhQEAjGEQELe+eOcC51K973dt/tWNcnv8zz93FtVp6re7tu33646myICMzOzgfaKqgMwM7OVgxOOmZmVwgnHzMxK4YRjZmalcMIxM7NSOOGYmVkpKk84kraVdI2khZIekXSKpFX62OdNkn6Tyz8naY6kCyRtXCg3WVLUeWwzsM/KzMyKVq3y5JKGAFcDdwNjga2AM0mJ8IRedl0XeBC4EHgE2AI4CRghaceIWFJTdhZwaGH/2c3Et/7668ewYcOaKTognnnmGdZaa63Kzt9It8YFjq0/ujUucGz9UXVcM2bMeDwiNqi7MSIqewBfAOYB69SsOxZYWLuuyWO9GwjgrTXrJgO39ze+ESNGRJWuu+66Ss/fSLfGFeHY+qNb44pwbP1RdVy9feZWfUttH+CqiJhfs24KMAjYtcVjPZF/rtaJwMzMrLOqTjjbkG55vSgi5pCucPqsZ5H0CkmrSdoaOA24Dbi1UGxbSfNzXc8NklpNZGZm1gGKCsdSk7QYOCYizimsnwtcGBHH97H/b4C98uIMYN+IeKxm+2eA50l1RBsARwMjgF0iopiYevaZAEwAGDp06IgpU6b056l1xIIFCxg8eHBl52+kW+MCx9Yf3RoXOLb+qDquMWPGzIiIkXU3NrrXVsYDWAz8Z531c4H/bmL/1wM7AR8mXSnNANbopfyapMYGlzYTn+tw6uvWuCIcW390a1wRjq0/qo6LLq7DmUdqcVY0JG/rVUTcGxG3RMRFpCudHYAP9VJ+IXAF8Nb+hWtmZv1VdcKZRaGuRtJmpCuRWXX3aCAiHgKeBLbsq2h+mJlZiapOOFcCe0lau2bdOOBZYHorB8oNB9Yj3TJrVGYQsB/p1puZmZWo0o6fwPnAUcAlkk4nXZ1MBM6KmqbSku4DpkfEYXn568AS4BbgKeCNpP4795OaVSNpXeBy4CLgPmB94LPAJsBBJTw3MzOrUWnCiYh5knYHzgMuIyWPs0lJp9aqQO1wN7cDR5Jak60BzAF+Dnw1Ip7JZZ4D/kkasWBDYBFwE7BrRNw+EM/HzMwaq/oKh4i4G9itjzLDCstTyFcyveyzCDig3fjMAIYd9+u2j3H08CWMb/M4s0/br+04zKpSdR2OmZmtJJxwzMysFE44ZmZWCiccMzMrhROOmZmVwgnHzMxK4YRjZmalcMIxM7NSOOGYmVkpnHDMzKwUTjhmZlYKJxwzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSmccMzMrBSVJxxJ20q6RtJCSY9IOkXSKn3s8yZJv8nln5M0R9IFkjauU3aspJmSFkm6W9K4gXs2ZmbWyKpVnlzSEOBq4G5gLLAVcCYpEZ7Qy67rAg8CFwKPAFsAJwEjJO0YEUvy8XcBfg58EzgK2Bf4saR5EfHbAXlSZmZWV6UJBzgcGAQcEBHzgamS1gEmSjojr1tORNwI3FizapqkucBvgTcDf8jrvwRcHxFH5eXrJL0JODGXNTOzklR9S20f4KpCYplCSkK7tnisJ/LP1QAkrQ6MAX5aKDcFGCVp3dbDNTOz/qo64WwDzKpdERFzgIV5W68kvULSapK2Bk4DbgNuzZu3Al5ZPD5wD+l5v6G90M3MrBWKiOpOLi0GjomIcwrr5wIXRsTxfez/G2CvvDgD2DciHsvbdgZuAHaIiD/V7PM64F5gr3r1OJImABMAhg4dOmLKlCn9fXptW7BgAYMHD67s/I10a1wwcLHNfPjpto8xdBA8+mx7xxi+aecvzFfGv2cndGtsVcc1ZsyYGRExst62qutw2nUk8Grg9aRGBldK2jkiFvX3gBExCZgEMHLkyBg9enQn4uyXadOmUeX5G+nWuGDgYht/3K/bPsbRw5dw5sz2/uVmHzK67TiKVsa/Zyd0a2zdGhdUn3DmkVqcFQ3J23oVEffmX2+R9DtSy7UPAd+r2b94/CE15zYzs5JUXYczi0JdjaTNgDVZvu6lVxHxEPAksGVedT+wuHj8vLwU+Gs/4jUzs36qOuFcCewlae2adeOAZ4HprRwoNxxYj3SVQ0Q8B1wHHFQoOg64KSLavylvZmZNq/qW2vmkDpmXSDqddHUyETirtqm0pPuA6RFxWF7+OrAEuAV4CngjcCzpqqa2lv9UUh+dc4BLSR0/9wX2HtinZWZmRZUmnIiYJ2l34DzgMlLyOJuUdGqtCtQOd3M7qcHABGANYA5pRIGvRsQzNce/QdL7gS8D/0Gu4/EoA7aiGNahxgztNoqYfdp+bcdhK76qr3CIiLuB3fooM6ywPIVlr2R62/dS0tWNmZlVqOo6HDMzW0k44ZiZWSmccMzMrBROOGZmVgonHDMzK4UTjpmZlcIJx8zMSuGEY2ZmpXDCMTOzUjjhmJlZKZxwzMysFE44ZmZWCiccMzMrRcujRUsaDrwN2Ig0NcCTpNkzb4wIT9tsZmZ1NZVwJG1Jmk/mEGAoaYrmp4DngFeRpoReKmk6cAHwk4hYOiARm5nZy1Kft9QkXQDcBWwPnALsAKwRERtExGsiYjCwIfBeYCZwBnCPpF0GLmwzM3u5aeYK51lgm4h4qFGBiHgcuBK4UtLngIOATTsTopmZrQj6vMKJiCN7SzZ1yi+NiJ9ExE+aKS9pW0nXSFoo6RFJp0hapY99dpT0fUn35f3+IukkSWsUyk2UFHUeezf7fMzMrDPammJa0nbAroCA6RExs8X9hwBXA3cDY4GtgDNJifCEXnYdl8ueDtwLvBk4Nf88sFD2aaCYYO5pJU4zM2tfvxOOpP8AvgJcA6wFfE3S0RHxzRYOczgwCDggIuYDUyWtA0yUdEZeV89p+TZej2mSFgHflrR54YpsSUTc3EJMZmY2AJppNLBmg03/BYyKiIMiYl/gCOCLLZ5/H+CqQmKZQkpCuzbaqZBsevwx/9ykxRjMzKwEzXT8/KukQ+qsF6l5dI/+NIPeBphVuyIi5gAL87ZWjMox3F9Y/ypJj0taLOmPkg7oR5xmZtYmRUTvBaR3AecAzwNHRcStef2nSc2kryH1w9kdODYizm365NJi4JiIOKewfi5wYUQc3+RxNgLuAK6IiPE16z9MarL9R2Bt4JPAvsCBEXFJg2NNACYADB06dMSUKVOafTodt2DBAgYPHlzZ+Rvp1rhg4GKb+fDTbR9j6CB49Nn2jjF803WXWe7WuDplZXyvtavquMaMGTMjIkbW29ZnwgGQJOAwUsX8VOC/IuLvkt7CS7e+ro+IP7USWCcSjqTVSA0PXgOM6G20g/w8bgQGRcT2fR175MiRcfvtt/dVbMBMmzaN0aNHV3b+Rro1Lhi42IYd9+u2j3H08CWcObOtdjrMPm2/ZZa7Na5OWRnfa+2qOi5JDRNOU2OpRXIBsDXwKHCnpC8CsyLif/KjpWSTzQPqfTUakrf1KieQC4E3Afv2NbROpOx6CfDmvppem5lZZ7U0eGdEzI+IY4CdSOOpzZL0/jbOP4tCXY2kzUi36GbV3WNZ55CaU4+NiGbKA0R+mJlZiZpqpSbpy5JuyZXuk4BFETGWVNdxkqTp+fZaq64E9pK0ds26caTRDab3EdcXSC3jPhwRNzRzsnxFdCDw54h4oR/xmplZPzVz4/a7wLakPjcLSUlmqqRtI2JqTjSfyusujYgJLZz/fOAo4BJJpwNbAhOBs2qbSku6j9Sx9LC8/CHgv4HJwMOS3l5zzPsj4p+53HTg56SrpbWAT5CuzvZvIUYzM+uAZhLOPsBBETEVQNLvgSdIPf3vy6NCnyfpR8BJrZw8IuZJ2h04D7iMNAL12aSkU4yzts5lz/xzfH7UOpSUiADuA/4T2JjUZPoPwH4RcWUrcZqZWfuaSTizgI9ImgEsIjUtfgaYW1soIp4EPtNqABFxN7BbH2WGFZbHs3yiqbffYa3GY2ZmA6OZhPNR0hXD46TK9tmkK55FAxeWmZmtaPpMOBHxF2CUpLWA1Tyrp5mZ9UfTvb0i4hnSrTQzM7OWNdMs+iOtdpKU9DpJ7+x/WGZmtqJppuPn54D7JZ3aW18bSetJOkTSZcCfSC3DzMzMgObqcHaQNA44EviipAWkCcweB54DXgVsAbyWNBzNRcDhEfHwgEVtZmYvO03V4eTpon8iaStgD+CtwEakzpSPAtcDvwemRcTiAYrVzMxexloaIjYi7mf5+WbMzMz61NLgnWZmZv3lhGNmZqVwwjEzs1I44ZiZWSlaSjiS3ivJScrMzFrWavK4FJgr6XRJbxyIgMzMbMXUasLZCvgOcDBwp6SbJH1C0jqdD83MzFYkLSWciJgdESdFxBbAu0kTnJ0N/F3SDyWNGYggzczs5a/f9TERcW1EfAR4AzADOAS4WtIDkj4rqaVOpWZmtmLrd8KRtKukycBfgO2A/yVN/XwxcDJwYScCNDOzFUOrrdQ2l3SipPuBa4HNgAnAxhFxZERcExHHkmYJHdvkMbeVdI2khZIekXRKX9MhSNpR0vcl3Zf3+4ukkyStUafszpJukbRI0oOSjmrlOZuZWWe0etvrAeAR0pTT34uIBxuUuwu4ta+DSRoCXA3cTUpQWwFnkhLhCb3sOi6XPR24F3gzcGr+eWDN8V8HXAVcDnwBeBtwlqSFEXFBX/GZmVnntJpw3gNcFRFLeysUEX8FmmlAcDgwCDggIuYDU3OLt4mSzsjr6jktIh6vWZ4maRHwbUmbR8RDef0xpAT54YhYAlwr6bXASZK+GxHRRIxmZtYBrdbh7EialmA5kjaWdGKLx9uHlMBqE8sUUhLatdFOhWTT44/55yaF41+Sk03t8V9DqncyM7OStJpwTiJ9WNezSd7eim2AWbUrImIOsDBva8UoYCl5+gRJa5HqmGYVyt1Tc24zMyuJWrmrJGkpsFNE3FZn21jguxGxfgvHWwwcExHnFNbPBS6MiOObPM5GwB3AFRExPq/bFJgLvC8iLq0puyqwGPhkREyqc6wJpIYQDB06dMSUKVOafTodt2DBAgYPHlzZ+Rvp1rhg4GKb+fDTbR9j6CB49Nn2jjF803WXWe7WuDplZXyvtavquMaMGTMjIkbW29ZnHY6kj5JanQEE8C1JxbqVNYDhwG/bCbQ/JK0G/BRYAHy23ePlJDQJYOTIkTF69Oh2D9lv06ZNo8rzN9KtccHAxTb+uF+3fYyjhy/hzJntdU+bfcjoZZa7Na5OWRnfa+3q1riguUYDC4En8u8CngaeLJR5HrgS+GaL558H1PtqNCRv65Ukkfr7vAnYOSJq93kq/ywef0jNuc3MrCR9JpyI+BnwMwBJ3wdO6aU5dKtmUahLkbQZsCbL173Ucw6pOfW7I6JYF/SMpL8Vj1+z3MzxzcysQ1odS+3QDiYbSFdFe0lau2bdOOBZYHpvO0r6AnAEqcnzDb0c/32FjqTjgL8Bd/Y7ajMza1nVc9ucDzwHXCJpj1xhPxE4q7apdB5R4Ls1yx8C/pt0O+1hSW+veWxQc/yvkVrV/VDSGEnHAp8kXaW5D46ZWYmaaTRwKzA+Iu6WdBup4UBDEfG2Zk8eEfMk7Q6cB1xGqnc5m5R0inHWXqXsmX+Oz49ah5JGQiAi7pO0N3AW6WrnH8DRHmXAzKx8zTQauIt0i6vn945eGUTE3cBufZQZVlgez/KJptG+N5CGtDEzswo102jg0Jrfxw9oNGZmtsKqug7HzMxWEs3U4fRZb1OrlTocMzNbeTRbh+MWXWZm1pZm6nDGlxCHmZmt4FyHY2Zmpai0H46Zma08Ku+HY2ZmKwf3wzEzs1K0PAlGnn9mPKn3/sbA34FbgB9ExPMdjc7MzFYYLTUakPRG4F7gf4HtgBfyz/8F7pO0bccjNDOzFUKrVziTSBOwvTMi5vSslPRa4HLS6M/v6lx4Zma2omg14YwEPlibbAAiYo6kk4AfdSwyW+kM69B0ye1Ouzz7tP3ajsPMltdqP5zZwBoNtq0BzGmwzczMVnKtJpzjgC9L2ql2paS3A6cC/9WpwMzMbMXSn8E71wFulPQY8BiwYX48ARwPXDoAcZqZ2ctcfwbvvGuAYjEzsxVY5YN35qbU5wKjSFNMXwCcHBEv9LLPasBXgLeTGjKsERGqU24y8NE6h3hjRMxqP3ozM2tWyx0/O0nSEOBq4G5gLLAVcCapbumEXnZdE/g4cCtwI71PUT0LOLSwbnb/IjYzs/6qNOEAhwODgAMiYj4wVdI6wERJZ+R1y4mIpyS9OiJC0hH0nnCeiYibOx+6mZm1ouXpCSSNk3S1pDmSHis+WjzcPsBVhcQyhZSEdu1tx4jwIKJmZi8jrQ5t8yHgB8B9wGuAX5FGGHgFMB84r8Xzb0O65fWi3Kl0Yd7WCdtKmi/pOUk3SOo1kZmZ2cBQKxcKkv4IXAycBiwGRkbEHyStDUwFLo6Ir7dwvMXAMRFxTmH9XODCiDi+iWMcAZzboNHAZ4DnSXVEGwBHAyOAXSLi1gbHmwBMABg6dOiIKVOmNPt0Om7BggUMHjy4svM3MlBxzXz46baPMXQQPPps3+V6M3zTdZdb162xdWtcndKt/wPQvbFVHdeYMWNmRMTIettarcN5PfD7iHhB0gukPjlExL8knQ6cDTSdcAZaRHyjdlnSFaRm3ccD+zfYZxJpzDhGjhwZo0ePHuAoG5s2bRpVnr+RgYqr3SFpIA1tc+bM9qomZx8yerl13Rpbt8bVKd36PwDdG1u3xgWt1+HMB1bPvz8MvLFmm4D1WjzePKDeV6MheVtHRcRC4ArgrZ0+tpmZ9a7VrzW3AW8GriLV35woaQnpttWJQKutwWZRqKuRtBmp2fNA9ZMJPGupmVnpWk04XwU2z7+fmH//FulK6Tbgky0e70rgGElrR8S/8rpxpCmtp7d4rD5JGgTsB8zo9LHNzKx3LSWc3J/l5vz7U8BYSasDqzfqM9OH84GjgEtyHdCWwETgrNrjSboPmB4Rh9Ws2wdYC9g+L78/b7otIh6StC6pBd1FpFZ16wOfBTYBDupHrGZm1oaOTTEtqeUppiNinqTdSc2pLyMNbXM2KekU41ylsO5bvHS1BfCz/PNQYDLwHPBP0ogFGwKLgJuAXSPi9lbiNDOz9rWUcPIU078hXSXMII0WvR3w78CXJO0dEXe3csxcvreRAoiIYc2sK2xfBBzQSixmZjZwPMW0mZmVotVm0SOBE+tNMQ2cBOzYqcDMzGzF4immzcysFK3eUjsOOFPSgxFxS8/KmimmP9/J4Mzs5WtYh0ZBaHc0hdmn7dd2HNYZnmLazMxK4SmmzcysFJVPMW1mZiuHfg0RK2kTYBTwatKttJsj4pFOBmZmZiuWVjt+rgKcC3yCZXv+vyBpEnBkRCztYHxmZraCaLVZ9MnAx0iNA4aRpoIelpc/xvJD0piZmQGt31L7d+CEwqyec4CvSQrSQJwndio4MzNbcbR6hbMhcEeDbXfk7WZmZstpNeH8FfhAg20fAP7SXjhmZraiavWW2peBKXmwzouBR0lXNQcBY2icjMzMbCXX6gRsP5X0FKnxwDeAVwKLSVMV7B0RUzsfopmZrQiaTjiSXkmadO3OiBgl6RWkWTQfd1NoMzPrSyt1OC8A1wLbAETE0oh4zMnGzMya0XTCyYnlXmCjgQvHzMxWVK22UvsicKKk4Z0KQNK2kq6RtFDSI5JOySMa9LbPapK+Jul3kp7NfYAalR0raaakRZLuljSuU7GbmVnzWm2ldgKwHvAnSQ+TWqkt82EfEW9r9mCShgBXA3cDY4GtgDNJifCEXnZdE/g4cCtwI7Bbg+PvAvwc+CapU+q+wI8lzYuI3zYbp5mZta/VhHMXcGcHz384aXicAyJiPjBV0jrAREln5HXLiYinJL06IkLSETRIOMCXgOsj4qi8fJ2kN5FGQ3DCMTMrUavNosd3+Pz7AFcVEssU4HRgV+CyXmJpeBsNQNLqpL5BRxU2TQG+L2ndiHi6X1GbmVnLmqrDkTRI0oGSjpb0IUlDO3T+bYBZtSsiYg6wMG9rx1akfkKzCuvvIT3vN7R5fDMza4H6uFBA0pakepZhNavnAwe3Ww8iaTFwTEScU1g/F7gwIo5v4hhHAOdGhArrdwZuAHaIiD/VrH8dqbXdXvXilzQBmAAwdOjQEVOmTGn9iXXIggULGDx4cGXnb2Sg4pr5cPsXnEMHwaPPtneM4Zuuu9y6bo2tW+OC7o6tE1a2/89mjRkzZkZEjKy3rZlbamcAS4F3kkYU2IJUCf/t/PsKJSImAZMARo4cGaNHj64slmnTplHl+RsZqLjGH/frto9x9PAlnDmzX/MKvmj2IaOXW9etsXVrXNDdsXXCyvb/2QnN3FIbRZqS4PcRsSgi7gE+CbxW0sZtnn8eUO/rx5C8rd1jU+f4QwrbzcysBM0knI2BBwrr7gdE+51AZ1Goq5G0GanZc7HupVX3k8Z5K9YFbUO6Yvtrm8c3M7MWNNvxs/eKnv67EthL0to168YBzwLT2zlwRDwHXEcaybrWOOAmt1AzMytXszdHr5K0pM76a4rrI6KVSdjOJzVbvkTS6cCWpGmqz6ptKi3pPmB6RBxWs24fYC1g+7z8/rzptoh4KP9+KjBN0jnApaSOn/sCe7cQo5mZdUAzCefkgTp5RMyTtDtwHqnPzVPA2aSkU2tVoDjczbeAzWuWf5Z/HgpMzse/ISeiLwP/ATwIfMijDJiZla/PhBMRA5Zw8vHvpvFIAT1lhjWzrsG+l5KubszMrEKtDt5pZmbWL044ZmZWCiccMzMrhROOmZmVwgnHzMxK4YRjZmalcMIxM7NSOOGYmVkpnHDMzKwUTjhmZlYKJxwzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSmccMzMrBSVJxxJ20q6RtJCSY9IOkVScTrpevutK+n7kuZJelrS/0lar1BmsqSo89hm4J6RmZnV0+cU0wNJ0hDgauBuYCywFXAmKRGe0MfuPwXeAHwcWAqcTppK+p2FcrOAQwvrZrcTt5mZta7ShAMcDgwCDoiI+cBUSesAEyWdkdctR9IoYE9g14i4Pq97GLhF0h4RcXVN8Wci4uaBfRpmZtaXqm+p7QNcVUgsU0hJaNc+9nu0J9kARMStwIN5m5mZdZmqE842pFteL4qIOcDCvK3p/bJ76uy3raT5kp6TdIOk3hKZmZkNEEVEdSeXFgPHRMQ5hfVzgQsj4vgG+00l3Srbv7D+ImDLiHhHXv4M8DypjmgD4GhgBLBLviKqd+wJwASAoUOHjpgyZUobz7A9CxYsYPDgwZWdv5GBimvmw0+3fYyhg+DRZ9s7xvBN111uXbfG1q1xQXfH1gkr2/9ns8aMGTMjIkbW21Z1Hc6Aiohv1C5LugK4Czge2L/BPpOASQAjR46M0aNHD3CUjU2bNo0qz9/IQMU1/rhft32Mo4cv4cyZ7b2tZx8yerl13Rpbt8YF3R1bJ6xs/5+dUPUttXlAva8fQ/K2ju4XEQuBK4C3thCjmZl1QNUJZxaFOhdJmwFrUr+OpuF+WaO6nVqRH2ZmVqKqE86VwF6S1q5ZNw54Fpjex34bSdqlZ4WkkcCWeVtdkgYB+wEz2gnazMxaV3XCOR94DrhE0h65wn4icFZtU2lJ90n6bs9yRNwE/Ba4UNIBkvYH/g+4oacPTh6J4HeSPilpd0njgOuATYD/LusJmplZUmmjgYiYJ2l34DzgMuAp4GxS0qm1KlAc7mZcLvs9UuK8HDiqZvtzwD9JIxZsCCwCbiJ1Fr29o0/EzMz6VHkrtYi4G9itjzLD6qx7ijRkTXHYmp7ti4ADOhCima1ghnWoBV07LfFmn7Zf2zG83FR9S83MzFYSTjhmZlYKJxwzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSmccMzMrBSVjzRg5eqGHtawcvayNlvZ+QrHzMxK4YRjZmalcMIxM7NSOOGYmVkpnHDMzKwUTjhmZlYKJxwzMyuFE46ZmZWi8o6fkrYFzgVGAU8BFwAnR8QLfey3LnAOsD8pcV4OHBURTxTKjQW+DLweeCAf+yedfh5mZu1a0TtmV3qFI2kIcDUQwFjgFOBo4OQmdv8pMBr4ODAe2BG4tHD8XYCfA9cB+wC/Bn4sac+OPAEzM2ta1Vc4hwODgAMiYj4wVdI6wERJZ+R1y5E0CtgT2DUirs/rHgZukbRHRFydi34JuD4ijsrL10l6E3Ai8NuBe1pmZlZUdR3OPsBVhcQyhZSEdu1jv0d7kg1ARNwKPJi3IWl1YAzpSqjWFGBUviVnZmYlqTrhbAPMql0REXOAhXlb0/tl99TstxXwyjrl7iE97zf0I14zM+snRUR1J5cWA8dExDmF9XOBCyPi+Ab7TQWeiYj9C+svAraMiHdI2hm4AdghIv5UU+Z1wL3AXhGx3G01SROACXlxa+Av/X6C7VsfeLzC8zfSrXGBY+uPbo0LHFt/VB3X5hGxQb0NVdfhdJ2ImARMqjoOAEm3R8TIquMo6ta4wLH1R7fGBY6tP7o1Lqj+lto8oF5dypC8rZ39en4Wyw0pbDczsxJUnXBmUairkbQZsCb162ga7pfV1u3cDyyuU24bYCnw137Ea2Zm/VR1wrkS2EvS2jXrxgHPAtP72G+j3M8GAEkjgS3zNiLiOVL/m4MK+44DboqIp9sPf8B1xa29Oro1LnBs/dGtcYFj649ujavyRgNDgLuBO4HTSQnjLOCciDihptx9wPSIOKxm3VWk0QM+T7piOR14LCLeWVNmF2AacB6pU+i+ufze9RoMmJnZwKn0Cici5gG7A6sAl5FGGDgbOKlQdNVcptY40lXQ94ALgRnA+wrHvwF4P7AHcBXwb8CHnGzMzMpX6RWOmZmtPKquwzEzs5WEE1I1XdsAABkPSURBVI6ZmZXCCcfMzErhkQa6gCSRGjzsB7wReDWp5d0/gJuByRFRSb+h3C9qX0DAzyLiCUmvIbX22wqYDUyKiJklxvRfwBVlnrNZkgYBq0bEv2rWbQAcAWxL+rv+Cfjmy6RpfiXy/8R7gbeSpi+5nfQ374pK5zyq/ePAbrlxUlUx7AasBvw6Ip7J77VPk1r8PkD633ykivjqcaOBiuU3yBXACFKCeR7YlPRPdiXpjbM1cGpEnFpybG8DpgJrAUuAJ4G9crwvAHcB2wEbAXtExO9Kimsp6fWZBfwImBIR95dx7r5IugK4NyI+k5dHkf6OS0ktKUX6Wz9P+rC6q6S4dgAGRcSNNev2Br7AS4nwz8DE2jIlxXYjcFhE3JOXh5CmDxkBLMjFBpO+fO1Vm8wHOK5P9bJ5EPA14BuksRmJiG+WERe8OCbkNcBmedWDpClbpgKvInV835rUp3FERMwtK7ZeRYQfFT6AH5PesMNr1m0C/Ab4eV7elfSP97GSY5tK6jz7KtLI2+cBc4FfAq/MZVYnfaBeV2JcS4GvkmZ5fY6U/G4DPgtsWvHf83FgbM3yzaQPhrVr1q1LatJ/VYlx3Qx8sWb5Y/l1vAb4InBC/lsvqY2/xL/n22qWv0v6crN3zbq9ScNRnV1yXC/kn/UetdteKPk1+ynpC8LrSHdEfpg/R27sea+RBvH8M/DtMmPrNe6qA1jZH6RptQ+ss35YfkNvnJePB/5ccmxPAPvULG+Y/7n2LJTbD3i8xLhe/IDK/2wT8gfnkvyYlte9uoK/50LgXTXLzxdfr5rX7JkS45pfGwdwH3BunXLnV/A+KyacfwL/Wafc54GHSozrF8DfgUPJd4Nqtr0qx/2usuIpnP8R4OCa5c1zPAcUyh0K/LWKGOs93Gigeq8gJZaiF0i3X3oGH72FaubwiTq/F+/DVnZfNiKejIhJEbE78BrSFOWrkT44/y6p/UniW3MnaeK/Ho+SkmLReqTkVJalheXNgYvrlLuYdCumSq8i1dkUzSDdvi1FRLwP+ChwDHBbnvLkxc1lxdHAENIt+B4P558PFco9QPq/6ApOONWbCnxZ0pY9K/I97P8hvaF6GgsMBsquZL4d+LykwZJeQbrKehj4D0mr5FhXBT5F+qCtVET8IyK+ERHvALYgjVixSclhnAYcJ+lj+bX5CvA1Se+WtJqk1XPdyVdZfjbagfQ74JCa5buAekPY78hLH15lOlDSp3K9yTyg3nwq65Ou1EoTaVSSN5NmCv61pCm50UzVHiN9aejxAvBt0hecWhsCpdR5NaXqS6yV/UH69nEXaWTr+0hjyz1LutVWezvrDOAnJcc2kvTPvzjH9ATwFlISfIA0HNGDpHqUMSXGtcwtmG57AB8nfTA+Ddyaf3+BdLvvhfz4BbBmiTENz3H8EHgbaSr2x0gJ8d2kCufTgEXUuZ1Vwt+z+PhenXLfBn5X4d91I9IwWguAM/PfsapbapfWe43qlDsXmFrVa1Z8uJVaF8hXCweTPszXICWeH0XEk5UGBuRvc+8hNaH/eUT8XdJGwLGkWy8PARdExB9KjOkk4DvRRc09iyStRxrv722kDyqRkvc9wOURMaOCmLYHvgXsRLolpLyp5/d5wCkR8Y2yY2uGpE8A90fEtRXHMYo05uPWwH5Rcqu+HMNQ0heWB/so9zlSndw15UTWOyccs5WMpDeSkk4xEd4YEYurjM1WbE44XUTSm0gTxNXOSjorSuqr0SpJr4iIYmV0ZSStQeqMuhS4r+oPz1yHsyU1HXkjYk6VMb3c5A6gRIUfVLkzryJiYc267ckdn6u4Wn25cqOBLpArmB8C7gB+RppAaVL+/Q5JsyUdWlFsB0i6VNIVkt6b142TNBtYLOmhfKujzJg+LOljNcurSjqN1HfjDlIDhiclHVdmXDXxjJD0K1Jl7T3A74GbgAclPSzpFElrVhFbN5K0Z2ESRiTtL+kPpPrD5yXdLmm/kuNaV9IvSHVf8yV9R9Iqkn4A/IH0/3mrpBskrV9mbM2SdKCkeq1gK+GEUzFJR5IqQy8HRpNalbwyPzYkdfq8HDhf0qdLju1gUjPZ9Un/+D/JyeWHpH4vR5E6mp0vaa8SQzue1OG0x+k5lq8C7yK9ZmcCJ0k6vsS4kLQn6TXZhHSf/1RSr/kXgImkCQYPBG7MrRHLjO09kq6RdI+kX0p6V50yO1XwAXUlaUinnhjeB1xCasBwXH4sBn6ZX9+ynAq8E/gcqaPsO0gtC3cjdUQdSqrf3CKXtT74llrFJD0AnB8RZ/RR7ljg8IjYsrdynSTpNmBGRByelw8hTXh3XkQcXVPu+8BmEbFHSXEtJLXgm56XHwO+UqzslvR54MiI2LzOYQYqthnAnRHx0cL6I0l9hLYk9RO6Ebg5InobPqWTcb2bNHrFzcAfgVHA9sA5wOd7bllJ2olUl1Oc8HAgY1sKvD0ibs3LfwAejoj3FspdAawVEbuWFNds0vvqO3l5B1JfoEMj4gc15T4BHB8RW5QRVz7n95osujkwusy/Z298hVO9jUlNZ/tyKyV2esu2ZtnOgZeTrryKnSkvIY3HVZanSVddPdYlDeFR9GfSVWKZtgUuqrP+IuC1wNYRsYj0Qf++OuUGyknAhRGxc0QcEREjgE8AnwQuyfVf3WI70lV/0STSYJ5l2YCX+sFBHjONNE5Zrfuo329oIH2U1JR9eB+P0r5sNcMJp3p/Bj6RO1bWlStOP0GqnyhTsOzU3j0DKT5VKLeA1Du8LL8idUhdLS9fDXywTrkPkvo4lekxUvP2oreQXs+ezrsP8dIoEmXYjkIijIjvkW4/vh24VlK9ERHKUnur5Wnqd1Z8hnI/sx4kvT493klq/PGOQrmdgbIbg9wLXBsRO/b2IN2O7BqenqB6R5Nuddwt6RLSCMg9H+jrklqtvY/UQXTvkmN7iPSN/SqAiHgh90G4p1BuK9KYU2X5Aqnn/J2SLiB1QD1d0nakcdREGl5mB9IQ92WaBJwqaTApET5P6r3/RdIApz19h7ak3A+pRaRRv5cRETPykC1XkW7zTSwxplpXSVqSf1+X9LebXiizDeW+z84HviFpOCkJHkx6730p/33/TLri+izl1+HczPKJr57a/laVc8KpWET8PjexPJY09MhmhSJ/I1Wqfi3KH4L/EgrjMEXELXXKfQAobU6QiHhS0ttJH+Kf46XbZqPy43nSkEHvjIjbyoorx/aVXCdxHHBiz2rSqOD/WVN0MfDfJYZ2B2l0gV8VN0TEAznpXAFMLjGmHifXWfdYnXUHkka0LkVEnJfvPHyQ1DDg2Ig4X9Jc0tBTPePhnQ98vay4snNJLeX6Mp1lx/arlBsNdJncXLbn9tRTtW3/u5Wk15JiLXWcq5rzD2PZToz3d0EfnFeSrvzWAB6o6rWpieeTpNZ9OzQawULSWqQhd/aICN9u70W+zb1+RPyz6lheTpxwzMysFL6l1iWUpnLehPTt/PE629cH9o2IC0sPro58D/sPwCFl37ZSl0/jrC6clrvb6WUyXXK+sqmd+noGKd7Sv7lLGkm6zSjSNPSzJL2FdIuy5312XkRcVXZsjfgKp2KSVie1Hjogr1pKGpH2c7UflhX1j9i3l81rAT8h1VXcCRARV5QUV1dO45xj6cppuZulNM7aQRFxSonn7MrpktWlU1/nWPYiNZZ5ktR6bwNgLKne9R5SX6sRpAYrB0bEpWXF1qsyhqT2o9fhw08ktUr7BGk6gKNIc1rcC7y+ptxOlD+NbVdOsUuXTuOcz9uV03K3EP+BFbzPunK6ZLp06ut83t+ThtZZJS8fn+P4bqHcD0kdjCt/b0V4iunKH6Rm0EcU1m0EXE+aandUXldFwrmdl6bY3bzweHP+hzy4Z12JcXXlNM75nN06Lfdrm3wcXsH7rCunS66TcLpi6ut8zqdJV8g9y0NyvLsVyu1JatBTWmy9PVyHU73NKHTojIh/SNqd9O3k6jykTJn9D3rsSLryOp3U7+XzkeffkNTTafEfEVGc1nag9UzjfH1e7pZpnHt047Tcs5s8p5os10kvl+mSu2Lq6+xZlu1X1fP7oEK5NUl9sLqCE071HgFez0sfngBEatb7AUnfIF06l95YINJXpEmSfgp8GZgp6dz8e5VOA/5P0t9Ir0vPNM5PkG6j9XT8LHsaZ3hpWu4bSMmudlruayN1nq1iWu5/AdcCF/RRbhdSn7AydfN0yQfmynnooqmvSbfUTpR0bz7310mzBf+XpOsj4l/5S+GxpITYHaq+xFrZH6TBMKf1UeYLlFxP0iCON5N68j8CfIZqp9jtummcc1zdOi33VNJQKH2Vq6IOpyunS6aLp74m1XfNrnmv30+6JdrzvzCTlJznAduXGVtvD1/hVO+bwDhJr44GHfIi4qtK8+W8u9zQlovjDmC0pA8AZ1DhkBkRcUGeq6RnGucn6YJpnCPi9jwUSnFa7nfx0rTcV1LytNykK+gJTZT7J4Wr7RJ8knTrpy8PkpJTKaL5zq+3k1psliYi7stDOe1MapxyTUQ8K2k06cvY1qRb8j+Kklr1NcPNoq1f8m2htYAFEdE1EzyZWfdywjEzs1J4vKSXiTy97XerjqOebo2tW+OC7o6tW0m6WtI1VcdR1K1xQffF5jqcl48xdO8XhG6NrVvjgi6NTdLVpDsfu1cdSx2iC18zujcu6LLYfEvNzF6Uvw2/IiK6Zkh7W3F0Teaz3klaI08D0HW6NbZujQu6N7aI2L1bk42kV3bja9atcUH3xeaE8/KxH6lZaDfq1ti6NS7o0tiq+oCS9GlJ90t6VtKfJX2kTrG3UvJr1q1xdXtsjTjhmK0kuvUDKvfrOpc0COuXSJ0YJ0u6WNIaZcbycoir22PrjRsNVEzStU0WrTekxoDq1ti6NS7o3thqPqB+TBq6/h2kD6ixwIcjosrxtj4PfD0iXhxSJ48l+H/AdZLeExFPOK6XTWwNudFAxSQtAf5CGgepN5sCO0W58+F0ZWzdGhd0b2ySbicNbVPvA+pB4D2RJoqrYt6lfwHvjYhphfXDSKMyrEKaBmCDMmPr1ri6Pbbe+AqnencBsyJiXG+FJL2fkofPoHtj69a4oHtj25r0rfhFEXGNpLeTPqBukrR3ifHUepo0AOYyImK2pHcAvwZuAk51XC/q5tgach1O9W4G3t5EuaD8scu6NbZujQu6N7aGH1Ck22uPkz6gdiwxph4zgP3rbYiIecDupPHK/qfMoOjeuKC7Y2vICad6ZwBHNlHuCmCLAY6lqFtj69a4oHtj6+YPqIuALSXVm9OIiHgW+DfS1ApzHBfQ3bE15Docs5WApIOAz5LqauqOSi5pFeBbwLsjouxEbSsBJxwzMyuFb6mZmVkpnHDMzKwUTjhmZlYKJxwzMyuFE471StJ4STMk/UvSPEl/lHTWAJ3rYEnjmyg3UVLUPB6R9HNJWzV5nsm5533lmn3OuWzP8763wfZ78/aJAxVDi8dd5nXu9HkkvULSEfk9+ayk+ZLukvQ/kvrVx0nJnyR9tMH2ybk3f71t58mT6vXKCccakvQFUjv+q4ADgH8Hfklq3z8QDgbGN1n2aWBUfnwe2B64RtJaTex7agvnGWitPGeARcAWkkbWrpS0IzAsbx/oGJpVfJ07fZ6fAF8GLiG9Jz9K6t/0juh/89uDgVcDP+rHvl8HDpH0un6ee4XnoW2sN0cA346I42vWXSbp5KoCqrEkIm7Ov98saQ7wO2Bf4GfFwrmPySoR8XxE3F9inJ32DPAH4AOkjpo9PgBcC4yoIqgeZb3OkvYB3g/sGxFX1mz6RX+vbrKjgB9GxOKac61KSp4fATYBPijpfuDkiHhxeKI8rMwNwH8AR7cRwwrLVzjWm1cB/yiurP322HPbRNL+kmZJWiTpBknbFvfLt1RmSnpO0t8kfSX/MyNpMnAgsGvNrbKJLcQ6I/8cVieuu0jf/Heq3VaI7V2SrpO0QNLTkqZJ2qFm+zslTZe0UNITkr4jae3eApI0StKvJP1d0jP5Vs0hta9dP5/zFODgng/W/PPgvL5jMeTX4OLC8UbnMtvVvpZ9vc6NziNpX0lLJW1ROM8Wef3YBq/BrvnncqNz9/fqJl+ZvAO4uLDpM8CxpFEYrgA+BnwPWK/OYX5OusrxZ2sdvsKx3vwBODJfPVzey3DnmwNnkebleBY4GbhK0ut7hr2XtCfpFsiFwDHAm0nfGtcDDs+/v5aU5D6Vjzu3hViH5Z//KKw7Azglr687z4uk0cBU4DrSbZlngJ1JIzr/UdLOwNXApaRv1esBpwFD8nIjmwO/B84nfRDvDHxf0tKI+DH9f86XkEYE2IV0VfdO0qjAlwBfKymGWsPo+3VudJ6/A4+QXveJNeXHA4+RBqGs55n882uSzoyIh1qMuZ7d83H/XFi/K2mk7TPyF6nf5zHo6rkRGAoMr3Mciwg//Kj7ICWFB0gDTS4ljYR8CrBOTZnJefs7atZtDiwBDq9ZdzNwXeH4xwIvAK/JyxcD05qIayJpsMlV8+MNpGQxH9i4ENf2dfafDNxes3wT6faUGpzvd3Vi3y0ff7smX0vlWL9N+vDqWd/Uc6593vn3XwL/m3//JnBp/v1xYGInYgCmARcX1o2ufd4tvs6NzvNlUpJSTZyzSfO9NHotNgLuyOcO4E7geGBwG+/3ScBtddZ/G/hbPudkYFgvx1g1v/c/0d84VuSHL/usoYi4A3gjqUL2m6QPgi8Bt0saXFP0sYi4sWa/h0i3uN4GL97XfyvL1638hHRbd1Q/wlsPWJwffwG2BMZFxN9ryjwcEX/q7SC5kcFOwA8if2IUtq+Z4/uppFV7HsAN+dwN60wkDVFqMfVQTawTSAmyXVOA90tanXSVtdzttBJi6NHn69yH75G+pIzOy2Py8vcb7RAR/wB2APYiXe29CvgKcKOk1QAkfTTfQvxTvo07K/8+Q9Ir6xx2I1LCLvoK6crnQdL/wufzVW+9uJYAT+VjWYETjvUqIp6LiMsi4oiI2Bb4OPB64LCaYo/V2fUxYOP8+/rAK4FHC2V6luuOeNuHp0lD6Y8EXkP61nlloUzxfPUMISXSv/eyfRVSwl1c83iO9Jw26+XYk4FxpNtce+Z4vwd0YgrgXwGDSR+GawGXVRBDj2Ze54Yi4gHS1dShedWhwK0RcVcf+70QEb+NiE+Rbtd9n3Qra1Te/oOI2J70ZWcJsHNEbB8RI6KmUUCNNUh/1+J55uTjvo90xb8LcIMadw94js6+visM1+FYSyLiu5LOALapWb1hnaIbkm7BQfrWuLhOuaH5Z93Ri/uwJCL66kvTTOXxPNLtwo0bbH8qH2ciqcK46JF6OynNK/8e4NMRcX7N+o58yYuIZyRdThoB+mcR8UyxTAdiWASsVlg3pF44TR6vNxcA31Fqin8ALbbyioilkn5LSlbFD/vXA/Oi7ymXn6TBlUlOUL9Rmqp7Immqh7MlnZMTUq1X0b/39ArPVzjWkKTlEomkDYB1WfZb7YZKswz2lHkt6VvlrZC+iZJusR1UONzBpA/7m/Ly85T8zTB/UN8C/HtPq686228Gto6I2+s86iYcYHXS/9eL35hzq7ZiH6Z2nvO3SFc25zfY3m4Mc1n2iwWkq6T+6u25XpK3TyHFXPcWIYCkoQ02/RuwkPT3rPUWmqvA/wt15iiq974Abss/X10ouwGwJvDXJs630vEVjvVmpqRfAr8l3SLbnNTJciHwg5pyjwMXSTqBl1qpPUa6ndPjJFLLte+TPkyGk1oufScielpFzQLGStqf9GH3SC8f6J10HKkV2pWSJpHu148iVXhfTmrccI2kpaSK73+RbuHsB3wxIpb7cImIpyXdBpwoaT4psR5HuhW4Tk3Rfj/nSPPZT+tle7sx/AI4TNLZpNZiY4B2pqFu+FwjYpGk/wM+Dfw4Ip7q5Tg/lfQv4KekxgUbAocAY0mV9cV930JqYNCX35Neqw0i4p81638k6Y/A9aTblyNIV5YPA/cUjjGSdMV3I7a8qlst+NG9D9I//29Jt40Wkf65fwRsU1NmMqmF1wGkb3XPkf5xl2u9RapLmEn6JjuXVP+was329Ukfck+Sb2M1iGsiubVWL7FPpqaFVF/bSE1frycl06dIrd62r9m+E/AbUku4Z4C7SU3B1+0lhtcB1+Tyc0iJa5nYm33OLTzvZVqptRsD8AVSC61/kWaZ/DeWb6XW1Ovc13MF9sjr9+jjOX4s/y3m5vfSk6SEOLpB+cuADzTxfl8NeAL4SGH9+/L5/kFK2vNJiX6HOsf4BoUWjX689PAEbNaW3KFvu4gY2VdZs97kusGDgS0jYmkHjzsH2Csiilcj9cp+A3hdROzXYPtkUqKcXWfbKsBDwHERcVFbQa+gfEvNzColaWtgW9KQMCd3ONkMIXWKbbZO5WvAXyW9IercKu3DQaRbyg3rn1Z2bjRgZlX7NulW7RWk4WM6JiLmRcSgSA1Xmik/l3TLrlGrxUtJt1zrEXBYpL44VodvqZmZWSl8hWNmZqVwwjEzs1I44ZiZWSmccMzMrBROOGZmVgonHDMzK4UTjpmZleL/AW46+1xJ9j1lAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot probability distribution\n",
    "x = uncertainty_model.values\n",
    "y = uncertainty_model.probabilities\n",
    "plt.bar(x, y, width=0.2)\n",
    "plt.xticks(x, size=15, rotation=90)\n",
    "plt.yticks(size=15)\n",
    "plt.grid()\n",
    "plt.xlabel('Spot Price at Maturity $S_T$ (\\$)', size=15)\n",
    "plt.ylabel('Probability ($\\%$)', size=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Payoff Function\n",
    "\n",
    "The payoff function equals zero as long as the spot price at maturity $S_T$ is less than the strike price $K$ and then increases linearly.\n",
    "The implementation uses a comparator, that flips an ancilla qubit from $\\big|0\\rangle$ to $\\big|1\\rangle$ if $S_T \\geq K$, and this ancilla is used to control the linear part of the payoff function.\n",
    "\n",
    "The linear part itself is then approximated as follows.\n",
    "We exploit the fact that $\\sin^2(y + \\pi/4) \\approx y + 1/2$ for small $|y|$.\n",
    "Thus, for a given approximation rescaling factor $c_\\text{approx} \\in [0, 1]$ and $x \\in [0, 1]$ we consider\n",
    "$$ \\sin^2( \\pi/2 * c_\\text{approx} * ( x - 1/2 ) + \\pi/4) \\approx \\pi/2 * c_\\text{approx} * ( x - 1/2 ) + 1/2 $$ for small $c_\\text{approx}$.\n",
    "\n",
    "We can easily construct an operator that acts as \n",
    "$$\\big|x\\rangle \\big|0\\rangle \\mapsto \\big|x\\rangle \\left( \\cos(a*x+b) \\big|0\\rangle + \\sin(a*x+b) \\big|1\\rangle \\right),$$\n",
    "using controlled Y-rotations.\n",
    "\n",
    "Eventually, we are interested in the probability of measuring $\\big|1\\rangle$ in the last qubit, which corresponds to\n",
    "$\\sin^2(a*x+b)$.\n",
    "Together with the approximation above, this allows to approximate the values of interest.\n",
    "The smaller we choose $c_\\text{approx}$, the better the approximation.\n",
    "However, since we are then estimating a property scaled by $c_\\text{approx}$, the number of evaluation qubits $m$ needs to be adjusted accordingly.\n",
    "\n",
    "For more details on the approximation, we refer to:\n",
    "<a href=\"https://www.nature.com/articles/s41534-019-0130-6\">Quantum Risk Analysis. Woerner, Egger. 2018.</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:10.720023Z",
     "start_time": "2020-07-13T23:35:10.711632Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">     ┌───────┐┌────┐\n",
       "q_0: ┤0      ├┤0   ├\n",
       "     │       ││    │\n",
       "q_1: ┤1 P(X) ├┤1   ├\n",
       "     │       ││    │\n",
       "q_2: ┤2      ├┤2   ├\n",
       "     └───────┘│    │\n",
       "q_3: ─────────┤3 F ├\n",
       "              │    │\n",
       "q_4: ─────────┤4   ├\n",
       "              │    │\n",
       "q_5: ─────────┤5   ├\n",
       "              │    │\n",
       "q_6: ─────────┤6   ├\n",
       "              └────┘</pre>"
      ],
      "text/plain": [
       "     ┌───────┐┌────┐\n",
       "q_0: ┤0      ├┤0   ├\n",
       "     │       ││    │\n",
       "q_1: ┤1 P(X) ├┤1   ├\n",
       "     │       ││    │\n",
       "q_2: ┤2      ├┤2   ├\n",
       "     └───────┘│    │\n",
       "q_3: ─────────┤3 F ├\n",
       "              │    │\n",
       "q_4: ─────────┤4   ├\n",
       "              │    │\n",
       "q_5: ─────────┤5   ├\n",
       "              │    │\n",
       "q_6: ─────────┤6   ├\n",
       "              └────┘"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# set the strike price (should be within the low and the high value of the uncertainty)\n",
    "strike_price = 1.896\n",
    "\n",
    "# set the approximation scaling for the payoff function\n",
    "c_approx = 0.25\n",
    "\n",
    "# setup piecewise linear objective fcuntion\n",
    "breakpoints = [low, strike_price]\n",
    "slopes = [0, 1]\n",
    "offsets = [0, 0]\n",
    "f_min = 0\n",
    "f_max = high - strike_price\n",
    "european_call_objective = LinearAmplitudeFunction(\n",
    "    num_uncertainty_qubits, \n",
    "    slopes,\n",
    "    offsets,\n",
    "    domain=(low, high),\n",
    "    image=(f_min, f_max),\n",
    "    breakpoints=breakpoints,\n",
    "    rescaling_factor=c_approx\n",
    ")\n",
    "\n",
    "# construct A operator for QAE for the payoff function by\n",
    "# composing the uncertainty model and the objective\n",
    "num_qubits = european_call_objective.num_qubits\n",
    "european_call = QuantumCircuit(num_qubits)\n",
    "european_call.append(uncertainty_model, range(num_uncertainty_qubits))\n",
    "european_call.append(european_call_objective, range(num_qubits))\n",
    "\n",
    "# draw the circuit\n",
    "european_call.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:11.080005Z",
     "start_time": "2020-07-13T23:35:10.832762Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAE+CAYAAAB1DJw3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd5xU1fnH8c9XUGmKKIoYdVdNjD0xGMvPBlgpagQVFaOIihp7jUoMoGLBrrGh2FEs2FCKgIAtFjCWiBoFAVEsKIiwoJTn98e5I7PDbJnd2blzd5/36zWv3Xvn3pnvDsM8c++55xyZGc4551wuVos7gHPOueTx4uGccy5nXjycc87lzIuHc865nHnxcM45lzMvHs4553LmxcMVNUn9JVna7StJwyVtEVOe5pKGSfo+ytMrWn+SpM8lLZM0sYJ922f8LanbskL+DVGWLaPXdp2M9b2iTC0KncklS+O4AzhXDT8CB0a/bw5cDoyXtK2ZLSpwllOBg4BjgS+BaZI2BO4A/gU8Acyr4jF6AtPTluPobLUl0A+4H5iftv4FYDegLIZMLkG8eLgkWGZmb0S/vyFpFvAK0JnwYV1IWwGfmNnw1ApJewCNgHvN7P1qPMb7ZvbfugpYG2b2HfBd3Dlc8fPTVi6JpkQ/S6PTSP+S9ImksujU0W2S1k5tLOnxbKeSotM230haPVpuLemB6JRUmaSJknZK234GcAKwY9opp/6EQgbwXvqprFxFz/dkxrrUqa7touXSaPkISXdJ+lHSbEkDJK2Wse8OkkZImi9poaS3JO0nqT0wItrs8+jxZkT7rHLaqqrXJfXaSLpO0jlRnnnR6b1yp8Vc/eHFwyVRafTza6AZ4Vt/X6ATcCnQkfJHJEOAvSRtllohScBxwMNmtjRa/QxwAHA+0IPw/2OCpN9G9x8KjAQ+Jpza2Q24Bzgtur9ntO6FKvI3ktQ47VaT/4eDgIXAYcDDwD+j31N/31bAa0Bb4JQo+9PAJsA70d8I0C3KfGglz1XV65JyBLAP0Af4O9AVuLIGf5tLAjPzm9+K9gb0B+YSTrE2JpyrnwAsANpm2b4xsDuhHWHTaN1qwCxgQNp2HaNttouWD4yW907bpjnhFM5daevuByZnPGf79Meq5G9JbZd5uyK6fyLwZGWPTSicBjyYsd27wLC05UeB2UDTCrJ0jR6nNGN9r2h9ixxflxnANKBx2rqbgK/jfg/5rW5u3ubhkmA9YGna8iygh5nNAZD0V+Bc4HeED7aULYFZZrZC0n3AsZL6W/hk60UoAqm2h52Bb81sUmpnM1sk6Xlgjzz/PUcSPmhTvqrBY7yYsTwV2DRtuSPhqGpxDR47XS6vywQzS79ybCqwgaTVbeXRnasnvHi4JPgR2JfwDfhr4KuoACDpUOBBwtVOlwA/EE7VPA00SXuM+wintDpIehvozspTN0T7fJvlub8B1s3nHwN8aLVvMJ+fsfwL5f/e9YA5tXwOyO11yZZJwJqUL/6uHvDi4ZJgmZlNruC+w4E3zexvqRWS9s7cyMxmSBpHOOLYjHAq69G0TeYAG2R5/DaEglQIS4A1Mta1quFjfU/44K+tYnhdXBHyBnOXdE2BnzPW9axg2yGEI46/Ac+YWfo35TcJp1j2Sq2Q1AzoAryav7iVmk24FDjd/jV8rPHAEZKaVHD/L9HPiu5PKYbXxRUhLx4u6cYSrqTqK2lfSTcQrvjJ5hnCt/s/EU5j/crMxgCvA49JOk5SV8KVVU2Ba+ssfXlPA7+TdGP0twxkZefIXA0AWgIvS+oRPd4FknpH938S/TxZ0i6Sts/2IEXyurgi5MXDJd1dwPXAWcBTQAlwdLYNzexnYBTwBTAuyyZ/IRSjmwiX+groaGaf5T921nwvENptDiMUkhLC31WTx/qE0KA9l3A58dPR486M7p9JaPPpRrikd0T2RwJifl1ccVLU7uhcvSepMeHD814zuzTuPM4lmTeYu3pP0hrAHwhHJOsRjlacc7XgxcM1BBsBbxEuOT3ZzGbHnMe5xPPTVs4553LmDebOOedy1iBOW7Vu3dpKS0trtO+iRYto3rx51RsWiSTlTVJWSFbeJGWFZOVNUlaoXd4pU6bMNbP1s94Z9+Bahbi1a9fOamrChAk13jcOScqbpKxmycqbpKxmycqbpKxmtctLxiCg6Tc/beWccy5nXjycc87lzIuHc865nHnxcM45lzMvHs4553LmxcM55+qjoUOhtJS9O3aE0tKwnEcNop+Hc841KEOHQp8+UFaGAGbODMsAPSua7iY3fuThnHP1Td++UFZWfl1ZWVifJ148nHOuvpk1K7f1NeDFwznn6puWLbOv33TTvD2FFw/nnKtPHnsM5s+HRo3Kr2/WDAYOzNvTePFwzrn64uWX4dhjYc894Z57oKQEk6CkBAYPzltjOfjVVs45Vz9MnQqHHAKbbw7PPAPrrgu9ejFp4kTat2+f96fzIw/nnEu6OXOgUydo0gRGjQqFo475kYdzziXZTz9Bly7w/ffhtFUN5y7KlRcP55xLqqVL4fDD4f33YcQI+NOfCvbUXjyccy6JzOCUU2DMmNA43qlTQZ/e2zyccy6JLr8c7r0XLr0UTjih4E/vxcM555LmvvugXz847jgYMCCWCF48nHMuSV58MQxyuN9+cPfdIMUSw4uHc84lxbvvQvfusO228OSTsPrqsUXx4uGcc0kwaxZ07gytWsELL8Daa8cax6+2cs65YjdvXriaqqwMXnsNfvObuBN58XDOuaL2889w6KHw6afhstxtt407EeDFwznniteKFdCrF0yaFGYH7NAh7kS/8jYP55wrVhdfDMOGwdVXw9FHx52mHC8ezjlXjG67DQYNglNPhQsvjDvNKrx4OOdcsXn2WTjzTDj4YLj11tj6clSm4MVD0jaSxksqk/SVpMskNarGfjtJelHSD9FtnKRdCpHZOecK5s034aijYKed4NFHV50RsEgUtHhIagWMAww4BLgMOA+otH+9pE2i/RoDf41ujYGxkkrqMrNzzhXMZ59B167Qtm0YJbdZs7gTVajQV1udAjQFupnZAsKH/9pAf0mDonXZdAHWAg41sx8BJL0OzAU6A3fUfXTnnKtD330X+nKYwejRsMEGcSeqVKFPW3UCxmQUiWGEgrJ3JfutDiwDFqWtWxitK76Tgc45l4uyMjjoIJg9G55/Hn73u7gTVanQxWMr4OP0FWY2CyiL7qvI8Gib6yVtIGkD4EZgHvBEHWV1zrm6t3x5uAz3rbfgkUdg113jTlQthS4erYD5WdbPi+7Lysy+AjoA3YFvols34AAz+64OcjrnXN0zg7POCldX3XJL6EmeEDKzwj2ZtBS4wMxuylg/G3jQzC6pYL+2wMvAVFa2b5wG7Aj8X3T0krlPH6APQJs2bdoNGzasRpkXLlxIixYtarRvHJKUN0lZIVl5k5QVkpU3n1k3GTaMLe66i1k9ejD9lFPy8piZapO3Q4cOU8xsp6x3mlnBbsC3QL8s6xcRikpF+90AzABWT1u3BjATuKWq523Xrp3V1IQJE2q8bxySlDdJWc2SlTdJWc2SlTdvWR991AzMevQwW748P4+ZRW3yApOtgs/VQp+2+piMto3oMtxmZLSFZNgK+NDMlqZWmNkvwIfAFnWQ0znn6s6kSWEWwL32gvvvh9WS11+70IlHAQdIWittXQ9gMTCpkv1mAttJWiO1QtKawHaEIxLnnEuGqVPhL3+BLbaAZ56BJk3iTlQjhS4edwI/A09J2jdql+gP3GBpl+9K+kzSkLT97gE2Ap6W1EVSV+AZoC0wuGDpnXOuNr76KvTlaNIERo0KEzslVEGLh5nNA/YBGgEjCD3LbwT6ZWzaONomtd8U4EBCR8GHgAcJp7r2M7P36j65c87V0k8/QZcu8MMPMHIklCR7cIyCz+dhZlOBjlVsU5pl3XhgfB3Fcs65urN0KRx2GHzwQZhCdscd405Uaz4ZlHPO1SUzOPlkePFFGDIEDjgg7kR5kbwmfuecS5LLLoP77oN+/aB377jT5I0XD+ecqyv33gv9+8Pxx4fiUY948XDOubowZgz06QP77w933VWUEzrVhhcP55zLt//8JzSQb789PPEErL563InyzouHc87l08yZ0LkzrLtuuLJq7bXjTlQn/Gor55zLl3nzQifAxYth3DjYaKO4E9UZLx7OOZcPP/8chh2ZNi20d2y7bdyJ6pQXD+ecq60VK8JAhy+/DI8+Cu3bx52oznmbh3PO1dZFF8Fjj8E118CRR8adpiC8eDjnXG3cdhtcey2cdhpccEHcaQrGi4dzztXUM8/AGWfAIYfAzTfXu74clfHi4ZxzNfHGG3DUUbDzzvDII9CoUdX71CNePJxzLleffQYHHQQbbwwjRkCzZnEnKjgvHs45Vx1Dh0JpKXt37Ahbbw1LloQJndZfP+5ksfDi4ZxzVRk6NIxTNXMmMoNly8LtzTfjThYbLx7OOVeVvn2hrKz8uiVLwvoGyouHc85VZdas3NY3AF48nHOuKi1bZl+/6aaFzVFEvHg451xlHn0U5s9f9VLcZs1g4MB4MhUBLx7OOVeRiROhVy/Ye+8w/3hJCSZBSQkMHgw9e8adMDY+MKJzzmXz4YdhlNzf/haefhpatYLjjmPSxIm0bwADH1bFjzyccy7TV1+FeTmaNYORI0PhcOX4kYdzzqVbsCDMBDhvHrzySjhF5VbhxcM551KWLg1zj//3v2EK2T/+Me5ERcuLh3POAZiFXuRjx8K998IBB8SdqKh5m4dzzgEMGAD33w/9+8Pxx8edpuh58XDOuSFDQvHo3Rv++c+40ySCFw/nXMM2ejScfHI4TXXnnQ1qQqfa8OLhnGu43nkHDj8cdtgBnngCVl897kSJUfDiIWkbSeMllUn6StJlkqo1BZekbpLelrRY0veSRktqXteZnXP10IwZ0KULrLtuuLJqrbXiTpQoBb3aSlIrYBwwFTgE2AK4nlDE/lHFvicC/wIGARcArYCO+BVjzrlczZsX+nIsWQLjx0PbtnEnSpxCf/CeAjQFupnZAmCspLWB/pIGRetWIak1cCNwhpndnXbX03We2DlXvyxZEoYdmTYtXJa7zTZxJ0qkQp+26gSMySgSwwgFZe9K9jsi+vlAXQVzzjUAK1bAccfByy/DAw/AXnvFnSixCl08tgI+Tl9hZrOAsui+iuwCfAKcIGm2pKWS3pT0f3UX1TlX7/z97/D443DttXDkkXGnSbRCF49WwPws6+dF91VkQ+D3hHaRvwMHAYuA0ZLa5Dukc64euvVWuO46OP10OO+8uNMknsyscE8mLQUuMLObMtbPBh40s0sq2O9FYD+gk5mNjtatDcwE/mVml2bZpw/QB6BNmzbthg0bVqPMCxcupEWLFjXaNw5JypukrJCsvEnKCnWft/Urr7Btv37M3X13Puzff9WJnXLQkF7bDh06TDGznbLeaWYFuwHfAv2yrF9EKCoV7fcYsAJokrF+HDC8qudt166d1dSECRNqvG8ckpQ3SVnNkpU3SVnN6jjv66+bNWlituuuZosW1frhGtJrC0y2Cj5XC33a6mMy2jYkbQI0I6MtJMNHgKJbud0JRcU551b1v//BQQfBxhvDc8+F+TlcXhS6eIwCDpCU3hunB7AYmFTJfs9HPzukVkhqCbQD3st3SOdcPfDtt2FCJwlGjYL11487Ub1S6OJxJ/Az8JSkfaN2if7ADZZ2+a6kzyQNSS2b2WTgWWCIpOMkdQGeA5YCtxXyD3DOJcCiReGIY84ceP75MJWsy6uCFg8zmwfsAzQCRgADCJ3/+mVs2jjaJt0xwDPADcCThMLRMXpM55wLli2Do46CyZNh2DDYZZe4E9VLBR/aw8ymEoYVqWyb0izrFgKnRjfnnFuVGZx5JowYAbfdBgcfHHeiestH1XXO1R+DBsEdd4TOgH/7W9xp6jUvHs65+uGRR+Cii8IpqyuvjDtNvefFwzmXfBMmQK9e0L493HcfrOYfbXWtyldY0kuStop+P1bSenUfyznnqum//4VDD4Utt4Snn4Y114w7UYNQnfK8J7BO9Pt9hDk4nHMufl9+GeblaN4cRo6Eddapeh+XF9W52uoL4HBJCwk9ujeLfs8quprKOefq1oIFoXDMmwevvAKbbhp3ogalOsXjKuB24GzAgEcq2E7R/TUfccw556pj6VI47DCYOjVMIfvHP8adqMGpsniY2d2SngN+B7wMnEaYRtY55wrPDE46KcwCeN99sP/+cSdqkKosHpKOBV4ws1clDQCeNbOv6j6ac85l0a9fmAVwwIBwhZWLRXUazNMbyf8JbFx3cZxzrhL33AOXXw4nnACXrjKNjyug6hSPecBG0e+pdg3nnCusUaPglFPgwANDL3JlztDgCqk6DebjgIckfRIt3y9pUUUbm9nOeUnmnHMpU6bA4YfDDjuEOchXXz3uRA1edYpHb8JghFsBfwI+B76ry1DOOferGTOgSxdo3TpcWbXWWlXu4upeda62KgOuB5C0L9DXzHwCJudc3fvhhzCh0y+/hCFI2raNO5GL5DQku5ltVldBnHOunCVL4JBDYPp0GDcOtt467kQuTc6jh0naXNIdkj6Q9GX083ZJm9dFQOdcA7RiBRx3HLz6Kjz0EOy5Z9yJXIacjjwktQMmAEsI84p/A7QBugM9JXUws3fyntI517BceGFoGL/uOjjiiLjTuCxynUnwOuA/QKeoLQQASc2AkdH9lc4S6JxzlbrlFrj+ejjjDDj33LjTuArketpqZ2BQeuGAXxvVrwN8smDnXM099RScfXYYYv3GG70vRxHLtXgsBiqaz2Ndwuks55zL3euvQ8+esOuuMHQoNPIxVotZrsXjBeBqSXukr4yWrwJG5CuYc64BGDoUSkvZu0OH0CjesiU89xw0bRp3MleFXIvHucB0YJKkOZLekzQHmEToPHhevgM65+qpoUOhTx+YORNBuMLqxx9hzJi4k7lqyLWfx/fAHpIOBP4MtAXmAG+a2Yt1kM85V1/17QtlZeXXLVkS1vfsGU8mV225XqrbyMyWm9loYHQdZXLONQSzZuW23hWVXE9bfSlpkCTv6umcqzmzMO94Nj6dbCLkWjzuBA4D/ivpTUl9JK1dB7mcc/XZNdfAwoXQOOPkR7NmMHBgPJlcTnIqHmbW38w2B/YDPgFuAOZIGhoNmuicc5UbOhQuvhiOPjpMI1tSgklQUgKDB3t7R0LkPLYVgJm9ZGbHAhsCZwC/B8ZImiGpv6SNKn8E51yD9NJLcPzx0KED3HsvHHMMzJjBpJdeCkOve+FIjBoVjzQ7AXsR5vqYB7wCnAh8JumYWj62c64++eCD0HN8yy1DT/I114w7kauFmoyqWyKpn6RpwHjC5bq9gY3M7K9ACXAXcG1ekzrnkmv2bOjcGVq0CNPJrrNO3IlcLeV6qe4EYE/gS+A+4D4zm5m+jZktl/QIcFbeUjrnkmvBgjAT4I8/wiuvwCabxJ3I5UGuRx7fAp2B0qjxfGYF270LZJ04StI2ksZLKpP0laTLJFV7EBtJq0maLMkkdc0xv3OukH75Bbp3h6lTYfhw+MMf4k7k8iTXHuY9qrndUmCVwiKpFTAOmAocAmxBmOJ2NeAf1YxxIrBxNbd1zsXFDE46KcwCeP/9sN9+cSdyeZTrfB4ASNoY2BJoknmfmY2sZNdTgKZANzNbAIyN+on0lzQoWlfZ87YCBgIXAffUJLtzrkD++U948EG4/PIwK6CrV3Jt81gLeBzYP7Uq+mlpm1V2CqoTMCajSAwDrgH2pupReS8HXiM01DvnitXdd8MVV8CJJ4axqly9k2ubx1XApoRGcwGHAu2BIYRRdXetYv+tgI/TV5jZLKAsuq9CknYgXNV1fo6ZnXOFNHIknHoqdOoEd9zhEzrVUzKzqrdKbSxNJ7RNPAYsBXYxs7ej+64HNjGzCicclrQUuMDMbspYPxt40MwuqWTfSYTRey+UVEooVgeZ2fMVbN8H6APQpk2bdsOGDav235lu4cKFtGjRokb7xiFJeZOUFZKVN66sLT75hB3PPpuyTTbh3ZtvZnk15+Xw17bu1CZvhw4dppjZTlnvNLNq34BFwJ7R7z8BB6Tdtw8wv4r9lwJnZ1k/G7iykv2OBL4G1o6WSwmnyrpWJ3e7du2spiZMmFDjfeOQpLxJymqWrLyxZJ0+3axNG7PSUrM5c3La1V/bulObvMBkq+BzNdfTVl8AraPfPwXSL5XdhaqnoZ0HtMyyvlV03yokrU7ocHgNsJqkdYDUYIzNo3YY51ycvv8+nKb65ZfQCXDDDeNO5OpYrldbjQX2BZ4GbgQekNQO+JkwTMn1Vez/MRltG5I2AZqR0RaSpjnh0twbolu6YcA04LfV/xOcc3m1ZAkcckgYm2rcONiq0uZLV09Uq3hIakroHPgtMF1SGzN7SNJCwhDtTYHTCcOSVGYUcIGktczsp2hdD2AxYSrbbBYCHTLWbQg8ClwCvFSdv8E5VwdWrIC//hVefx0eewz22CPuRK5AqiwekjYndOwrTVu9QNIRZvY04Sikuu4EzgSeknQNsDnQH7jB0i7flfQZMMnMTjCzZcDEjEypLB+Y2Zs5PL9zLp/OPx+efBKuvx4OPzzuNK6AqtPmMQhYQbg8txmwLfAfqj7KWIWZzSM0rDci9OkYQDj91S9j08ZU3l/EORe3m2+GG2+Es86Cc86JO40rsOqcttoNOM/MXouWP5J0cvSzrZnNyeUJzWwq0LGKbUqruH8GKzsoOucKbfjwUDC6dQtHHd6Xo8GpzpFHW2B6xrpphA9vv6TCuYbmtdfCpE277QYPPwyN/CRBQ1TdS3Wr35PQOVd/ffIJHHxwmDL22Wehmp0AXf1T3Ut1x0halmX9+Mz1ZrZB7WM554rON9+EvhyNG4e+HK1bV72Pq7eqUzwG1HkK51xxW7QIunYNBWTiRNh887gTuZhVWTzMzIuHcw3ZsmXQowe88044VfXnP8edyBWBGs3n4ZxrIMzgtNPghRfgzjvD0Ydz5D4ku3OuIbnqKhg8GC6+GE4+Oe40roh48XDOZffww2Eip2OOgYED407jiowXD+fcqsaPh969oWNHGDLEOwG6VXjxcM6V98EHoef4738PTz0Fa6wRdyJXhLx4OOdWmj079OVYa60wnWzLbNPvOOdXWznnUn78ETp3hp9+gldegU02iTuRK2JePJxzYQbA7t3ho49g9GjYYYe4E7ki58XDuYbODE48MTSSP/AA7LNP3IlcAnibh3MN3aWXwkMPwRVXwLHHxp3GJYQXD+cassGDQx+Ok06CSy6JO41LEC8ezjVUL7wAp54aGslvv937cricePFwriGaPBmOOAJ23BEeeywMs+5cDrx4ONfQTJ8OXbrABhvA889DixZxJ3IJ5MXDuYbk++9DJ8ClS8OEThv6TNKuZvxY1bmGYvHiMIXszJkwbhxstVXciVyCefFwriFYvhz++lf497/h8cdhjz3iTuQSzouHcw3B+efD8OFwww1w2GFxp3H1gLd5OFff3XRTuJ19NpxzTtxpXD3hxcO5+uzJJ+Hcc8O4VddfH3caV4/4aSvn6pOhQ6FvX/aeNStcijt3Luy2Wxh+ZDX/rujyx4uHc/XF0KHQpw+UlSGAb74JvcaPPRaaNo07natn/KuIc/VF375QVlZ+nRlcdVU8eVy95sXDufpi1qzc1jtXC148nKsvKpr5b9NNC5vDNQgFLx6StpE0XlKZpK8kXSapURX7/FnSfZI+i/b7RFI/SU0Kldu5omYGW2yx6vpmzcKQ687lWUGLh6RWwDjAgEOAy4DzgAFV7NoD2AK4BugM3AacCwyts7DOJcmVV8KECWH4kZISTIKSkjBfR8+ecadz9VChr7Y6BWgKdDOzBcBYSWsD/SUNitZlc7WZzU1bnihpCXCXpBIzm1nHuZ0rXg89BP/4Rxh+5IEHQGLSxIm0b98+7mSuHiv0aatOwJiMIjGMUFD2rminjMKR8p/o50b5i+dcwowbB717h3nH77nHJ3RyBVPo4rEV8HH6CjObBZRF9+ViN2AFMC0/0ZxLmPffh27dYOutw7hVa6wRdyLXgMjMCvdk0lLgAjO7KWP9bOBBM6vWJMqSNgTeB0aaWa8KtukD9AFo06ZNu2HDhtUo88KFC2mRoMlykpQ3SVmhuPKu+e23/Om00wB45/bb+Xn99cvdX0xZqyNJeZOUFWqXt0OHDlPMbKesd5pZwW7AUuDsLOtnA1dW8zHWAF4GpgOtqrNPu3btrKYmTJhQ433jkKS8ScpqVkR5580z2247s7XXNnv//aybFE3WakpS3iRlNatdXmCyVfC5WugG83lAyyzrW0X3VUqSgAeBbYHdzazKfZyrV375JZyq+uSTMBPg9tvHncg1UIUuHh+T0bYhaROgGRltIRW4iXCJ735mVp3tnas/zELj+IQJ8OCDoZHcuZgUusF8FHCApLXS1vUAFgOTKttR0sXA6cAxZvZq3UV0rkj17RsGPxw4MFyW61yMCl087gR+Bp6StG/UqN0fuMHSLt+NepIPSVs+GriScMrqS0m7pt3KtxQ6Vx/ddVcY4LBPH7j44rjTOFfY01ZmNk/SPsC/gBHAfOBGQgHJzJU+ZMn+0c9e0S3d8cD9+U3qXBF5/nn429+gSxe47Tbvy+GKQsHn8zCzqUDHKrYpzVjuxapFw7n67+23oUcP+NOf4LHHoLFPweOKg4+q61yxmj4dunaFNm3C0Ufz5nEncu5XXjycK0Zz58KBB8KyZeGS3DZt4k7kXDl+DOxcsVm8OIyOO2sWjB8Pv/993ImcW4UXD+eKyfLlcMwx8MYb8MQTsPvucSdyLisvHs4VCzM491x46im46Sbo3j3uRM5VyNs8nCsWN94It9wC55wDZ50VdxrnKuXFw7li8MQTcN55cNhhcN11cadxrkpePJyL2yuvhOFG9tgjzAq4mv+3dMXP36XOxemjj+CQQ6C0FJ59Fpo0iTuRc9XixcO5uHz9NXTqFGYAHDUK1l037kTOVZtfbeVcHBYuDGNVzZ0LEyfCZpvFnci5nHjxcK7Qli2DI46A996D556DnbLP8ulcMfPi4QaGpP8AABm8SURBVFwhmcGpp4bTVIMHQ+fOcSdyrka8zcO5Qho4EO65B/7xDzjppLjTOFdjXjycK5QHHoBLL4Vjj4XLLos7jXO14sXDuUIYOxZOPBH23RfuvtsndHKJ58XDubr23nthnKpttoEnnwyX5jqXcF48nKtLX3wRGsVbtoSRI8NP5+oBv9rKuboyf37oBLhwIbz6KvzmN3Enci5vvHg4Vxd+/hkOPRT+9z8YPRq23z7uRM7llRcP5/JtxQro3Tv0HH/4YejYMe5EzuWdt3k4l299+8Ijj8CVV0LPnnGnca5OePFwLp/uuAOuvhpOOQUuuijuNM7VGS8ezuXLc8/B6adD165w663el8PVa148nMuHt96CI4+Edu1g2DBo7M2Jrn7z4uFcbU2bFo422raFESOgefO4EzlX57x4OFcbc+eGvhwrVoSRctu0iTuRcwXhx9bO1dTixXDwwaEX+fjxsOWWcSdyrmC8eDhXE8uXh8tw33gjjFf1f/8XdyLnCsqLh3O5MoNzzoGnn4abb4Zu3eJO5FzBFbzNQ9I2ksZLKpP0laTLJDWqxn4tJd0naZ6kHyUNlbReITI7V84NN4RLcc89F848M+40zsWioEcekloB44CpwCHAFsD1hCL2jyp2fxzYEjgRWAFcAzwD7FlXeZ1bxWOPwfnnhznIr7027jTOxabQRx6nAE2BbmY21szuBAYA50pau6KdJO0G7A8cZ2bDzexp4BhgD0n71knSoUOhtJS9O3aE0tKwXMySlDdJWWFl3g4dQl+O3/8+zAq4ml+s6BquQr/7OwFjzGxB2rphhIKydxX7fWNmL6dWmNlbwOfRffk1dCj06QMzZyIzmDkzLBfrh1yS8iYpK5TPm1o3axYMHx5nKudiV+gG862Al9JXmNksSWXRfSMq2e/jLOs/iu7Lr759oays/LqyMujVKwx2V2z+9z9Ytqz8umLNW1nWyy8Py2blf2ZbV92ftd3nu+9CH450ixeH94gPeugasEIXj1bA/Czr50X31WS/zbPtIKkP0AegTZs2TJw4sdoh9541i2yjEtmyZXy3/vrVfpxCWX/q1MTkrTRr27a/jgdl2caFSq2raJuM+6uzTVXP03bEiOx5Z81iUg7vqUJbuHBhTu/5uCUpb5KyQh3mNbOC3YClwNlZ1s8Grqxkv7HAM1nWPwy8XtXztmvXznJSUmIWvnuWv5WU5PY4hZKkvEnKapa8vJEJEybEHSEnScqbpKxmtcsLTLYKPlcL3eYxD8g2iXOr6L5871czAwdCs2bl1zVrFtYXoyTlTVJWSF5e5wqk0MXjYzLaKCRtAjQje5tGhftFKmoLqZ2ePWHwYCgpCac1SkrCcrGe405S3iRlheTlda5ACl08RgEHSForbV0PYDEwqYr9NpS0R2qFpJ0I7R2j6iIoPXvCjBlMeuklmDGj+D8skpQ3SVkheXmdK4BCF487gZ+BpyTtGzVq9wdusLTLdyV9JmlIatnM/g28CDwoqZukvwBDgVfNbFxB/wLnnHOFLR5mNg/YB2hEuCx3AHAj0C9j08bRNul6EI5O7gUeBKYAh9ZlXuecc9kVfGBEM5sKdKxim9Is6+YDx0c355xzMfLxFZxzzuXMi4dzzrmcydKHZKinJH0HzKzh7q2BuXmMU9eSlDdJWSFZeZOUFZKVN0lZoXZ5S8ws6zAVDaJ41IakyWa2U9w5qitJeZOUFZKVN0lZIVl5k5QV6i6vn7ZyzjmXMy8ezjnncubFo2qD4w6QoyTlTVJWSFbeJGWFZOVNUlaoo7ze5uGccy5nfuThnHMuZ148nHPO5cyLh3POuZx58XDOOZczLx7OOedyVvBRdV1+RDMwdgYEPGFm30vaGDgf2AKYAQw2sw/iSwmS/g6MjDtHdUlqCjQ2s5/S1q0PnA5sA6wA3gVuN7Mf40npXPz8Ut2IJBHmB+kCbA2sS/ig+Bp4A7jfzP4XX8KVJO0MjAWaA8uAH4ADgJHAcuBDYDtgQ2BfM3slpqhIWgEYYbrgR4BhZjYtrjxVkTQS+NTMzoqWdyPMVrmCMIeMgHbAL0BHM/swxqw7Ak3N7PW0dQcCF7Oy0L0H9E/fplhE/+cOAv5EeI9MJnzRKOoPJUlrE8aK6mhmr8adB37N1BFYA3jBzBZFX3pOI8y4Op3wZfKrvD1nkf87FUT0Io8kfCh8Tfhg+A3hDT2K8OL/HrjczC6PK2eKpLGEo8ZDgUWECbX+QvhwO8zMlkpaE3gGaGJmHWLMugK4Btge2I+Q+x1CIXnczL6MK1s2kuYCJ5jZs9HyG4TX+C+poxFJLYHngCVmdkCMWd8ARpjZwGi5N3APMAF4iVDo9gH2BLqn/qaYsr5OeF0/ipZbEWYHbQcsjDZrQfiidkD6kV8cJP2tkrubAtcCNwOfApjZ7YXIlY2k3wLjgU2iVZ8D+xO+YK4DTCN8fi0G2pnZ7Lw8sZk1+BvwKOFNsH3auo2A0cDwaHlvwpu8dxHk/R7olLa8AeFb5v4Z23UB5sacdQWwc/T7ukCf6I2+LLpNjNatG/frGmUsA/ZKW/4l83VNe20XxZx1QXo24DPg1izb3Qm8Vyzvg2h5COGI+cC0dQcC84Abi+B9sIJwFL+iglv6fctjzvo44Qjzt9H/sYeiz7PXgbWibVpH29yVr+f1BvOgE3CRpZ2Xt3B4dwrwF0ltzWwScCVwVkwZM1mW3zMPI4vqsNLMfjCzwWa2D7AxcB7hMPtOYI6kF2INGPwXSD9S+4bwHzLTeoRCE6cVGcslwJNZtnuS8M2zmBwMXGZmo1Mrot8HAt1iS7XSs8C3wAlAIzNbLXUjvB8EtI/WZU6ZXWh7AAPN7DMz+wH4B6Hd8zqLjuDMbC5wE+Xf27XixSNYjfBNItNywpukZbT8JrBloUJVYjJwvqQWklYDLgG+BE6V1AhAUmPgb4QPw6JjZl+b2c1m9n/AZoR57DeKORbA1cBFknpHr+FA4FpJ+0laQ9KaUbvCVYRvfHF6BeiZtvwhkG3o7T8T3h/FZB3C+zjTFEJbXazM7FDgOOAC4G1Ju6ffHU+qCrUinG5PSf1bZ85hNJ3wpS0v/GqrYCxwhaT3zWw6/HpO9hbCP0qqobwFUAxX2PQlZJ5HOPVTRmgsewL4VFKqwXwjwqmAomZmMwkf2lcXQZanJJ1B+JZ2I/AJ4cvDaMKHhqJNnyN8sMTpEuC16AvErYSG8gckrUs4HSjC++Js4KK4QqbpLilV3OYB2SYZak04HRc7M3tR0g6E1+8FSaMJVzPG2h6TxbeEo86U5cBdhKPmdBuQx+zeYA5El7iOIRxVzCSc594M+Bk4ysxGRdsNIsys1SOurClR5q6ELwDDzWyOpA2BCwmnKGYC95jZOzHGRFI/4G7L41UehSBpPaAHsDPhm7AIH3gfAc+b2ZQY4/1K0h+BO4BdKF/cUr/PI5weujmehEF04USm+82sd8Z2dwHbmNmehUlWPdH/rUGEU2p3EQpKBzN7OdZggKRngB8yX8ss290KbGVm++Xleb14BNHpniOAPwBNCI2Pj0TnEJ0rapK2JhSQzEL3upktjTNbLiSdBEwzs5fizpJNdOn2jYQvaF2sCC6BltQGaGZmn1ex3bmECyfG5+V5vXjUP5JWM7Ns3/SKhqQmhEa9FcBnxfgBF7V5bE5anx8zmxVvKueKgzeYZ5C0raTukk6Mbt0lbRt3rkySukl6RtJISQdF63pImgEslTQz+hYXK0nHRP0PUsuNJV1NuEzzfUKD/g+SiuGcPACS2kl6jnB++CPgNeDfwOeSvpR0maRmsYasRxSJO0c2kppm/ltL+mP0udAurlxFIc7rk4vpBvQmtBNku7Z7OWG4j+PjzhllPSLK9SrhksIy4CRCW80QQq/SR6PcB8ScdSpwatry9VHeS4HdCZcZ9id0YLqkCF7b/QltXZMJl2b3J3QU/SXKfB7hqqZ3gVZFkLcrod/MR9F7Ya8s2+xC/H0R9ifqc5C27i+EDqPLgKXRa94l7tc0ytYSeDrKtQy4G2gEPJDxufAq0DruvNX8m7rn830Q+x9UDDfgjOhNchuhN27r6I3SKPp9D+Bf0YfKaUWQ923gzrTlnlG26zO2uw8YF3PWMmDvtOVvgbOybHc+MLMIXtspwAMVvEdmEI7Wm0QferfHnHW/6APstej9OSVavp7olHS0XTEUj+WU7yR4aPQB/Hr0b38e4ehuGVk6ZcaQ9xbCECRnAMdGXxiGA19EhXB9Qv+wL4E74s5bzb8pr8XD2zwASdMJH8aDqtjuQuAUM9u8MMkqzLEA6GZm46LlloQG0n0traExOp11l5nF1n9C0hzgdDMbHi3/TDgampix3X7Ac2bWtPApy+VYDBxsZmMz1rci9Ozf1sw+knQscI2ZtY0jZ5TpVcI4XMenretN+OAbS7hScImkXQgN57F1ZouuttrVzN6Klt8BvjSzgzK2Gwk0N7O9Y4iZnmMGoePd3dHyjoTifLyZPZC23UmEI+bNYgkaMtxbzU1LCB0b8/I+8DaPoC3wVjW2e4si6MBEuAwz/Q2QGhtofsZ2CwmdseL0HKFD4xrR8jjgqCzbHUX4dhe3bwlX3GX6A+F1T/XzmcnKzqNx2Q54OH2Fmd1LGEpnV+ClqM9HMdqOcMlrpsGEgRLjtj4r+3dBNIYVYZyodJ+Rvb9KIR1HOBravopbSUUPUBPeSTB4DzhJ0stWwVVKUYPeSYRG3rjNJIyaOgbAzJZHlxB+lLHdFsCcAmfLdDGhJ/R/Jd0DjACukbQdKzuydQB2JIywGrfBwOWSWhAK3S+EHtp9gQm2sr/K5kDcV14tIYysXI6ZTYl6RI8hnBbqX+BcFUk/zfEj2TusLaI4vtR+TijCk6LlPQmn2f6P0M6Rsjvxvw8+Bd4ys2Mr20jSYcBj+XpSLx7BeYQexFMlPUUYPjz1Lb4lsBXhHO3GFEeP7afIGGbAzN7Mst2RlH+jF5yZ/SBpV8KH77mEXq4Au0W3XwinWPY0s7fjSbmSmQ2MTrFcBPwztZpwAcLZaZsuJTSox+l9wnn35zLvMLPpUQEZCdxf4FwVGSNpWfR7S8IXhkkZ22xF/F94IIy3drOk7QmF7gjCF6FLoy8W7xGOkM4B4h5p+w1CUatKeifSWvM2j4ikLQi9sw9k5dDGKV8Qrri51op4LopMkjYF5ptZUQz3ACCplPId2aZZcfbxWJ1w5NYEmF5Mr2GKpJMJQ5TsaBV0ZpXUnHDV0L4WBvWLRTTSQKZPzeyRjO0mRuuL4TLzMwmnU1cnjNZwp6SjCG1KqYExBwN/j/M9HF0yvLuZ3VLFdq0JbXaZBbtmz+vFY1XRdd2ptoL5Zhb36KnOuSIRncJubWbfxZ0lTl486pnokPodoGcxnAZSAqd1VUKm+HUuTl480kQfGhsRTqXMzXJ/a6CzmT1Y8HDlc3Su5O7mhEaxi4iGYzezkYXIlY0SNK0rJGuK3+qKxr063MwuizlHwadKzafoiCN92twphL8j9g9RhdGKuxP+P91vZh9L+gMwgJVfeP5lZmPy9qRxd1wphhuwJmE48+XRbSmhp3bLjO1i72wV5UjSLGdzgUPSlt8g9IheK21dS0LD6ZgieG3HEqZxXYdwrvtfwGxC7+3V094vowhXX8X+/q3G35TXzmE1zPBbwlWCqfflNMKH2nRCgX6bMBT7N8DGRfCavQ5snbbcKsq4Isq5gJWdHNeKK2eU7QDCl6+vo9d1AeEKxnlRvtui/3fLCdMp5+d54/5HKoYb4aqa+YRLcXcCzozexJ8Cv0vbrliKx2TCFSnHE67dTr/tEL2pj0itizlrYqZ1jXIkaYrfTat5OyXu9y0xTZVai7yJmTaXMMLAE4QZDyFcRDEPGJKx3UPAG3l73rj/kYrhRrg09/SMdRsCLwPfAbtF64qleIgw7/e3hCETNku7r2X0xl9ljKOYsr4F9Etb/gI4Mst2xwLfFUHe7zM+INaPXs/9MrbrXATFI3WUWdWtGI5AvwKOSFsuiXJ1y9jueOB/RfA+yCwe3wFnZ9ku9mF1CJcS75u23CrK3zFju/0JFwDl5Xm9n0ewCRmd/8zsa0n7EKr1OEk9KY7rz7HwThgs6XHgCuCDaKKXK+JNltXVwFBJXwAPsnJa1+8Jp6pSnQSLYVpXWDnF76uEo6b0KX5fstAhs1im+P0JeAm4p4rt9iBchh6nWKZKzaNinjZ3MeU7i6Z+zxzqpxmhY2leePEIvgJ+RzjS+JWFa7ePlHQz4bAw1obyTGY2Hzhd0mDCteefAtdQRHMsW7KmdYVkTfH7FqFd7oXKNormTolbLFOl1lJSps19DfinpE+jLNcRRrP+ezRqxk/R+HcXEopdXvjVVvw6sNjmZta+km0uJnxrNotxgLnKSDqSMFXmxoQB0GKfIjNFCZnWFRI1xe+lQB8zy+zUmrndXsAAM+tQmGRZM8QyVWpNKUHT5kr6LWEondT7YAbhaP5JQo/9mUAp4ctQBzN7Ny/P68Xj18vcegBXWSXTzko6mnDu+/iKtolbdEqlObDQzJbHncc5iG+q1LqmIpk2N+rftTvhCsHxZrY46ux8Iiu/8DxiZrPz9pxePJxzzuWqGEavdHVE0t2ShsSdozqSlBWSl9e5fPMG8xxIuhtYzcxOiDtLNXUgOV8QkpQVEpRX0jjCWYZ94s5SlSRlhWTlzXdWLx65ScwHBoCZ/TbuDNWVpKyQuLwiOe/bJGWFZOXNa1Zv86jHoks0NzCzuCerqVKSskLy8jqXb0mpmEVBUpNojoyk6EKYES0JkpQVEpRX0upJed8mKSskK2++s3rxyE1iPjBcwyDpNEnTJC2W9J6kv2bZ7E8Uwfs2SVkhWXnjyOptHgkkqbrXlGfrEVtQScoKycobdQq9lTBF7n8IU5HeL+kQ4Bgzy9tQFLWVpKyQrLxxZfXiQbI+MCJ7EYb5mFrFdsUwLEWSskKy8p4PXGdmv45bFY3HNhSYIKmrmX0fW7rykpQVkpU3lqzeYA5IWkb1PjB+A+wS9/Akkt4DPjazHlVsdxjwWJx5k5Q1ypGYvJJ+Ag4ys4kZ60sJ8400Ioy/tT7wumetviTljSurt3kEHwL/NbPDK7sBN8QdNPIGsGs1tksfeDAuScoKycr7I2FgvnLMbAbh1MVc4N/AnwsbK6skZYVk5Y0lqx958OvgZgeaWUkV23UnzGkda9GVtAWwrZk9V8V2TQmXk2YOe10wScoa5UhMXknPAj+Z2TEV3N+UMDheJ2Ie0DNJWaM8ickbV1YvHiTrA8O5FEmHA+cAXSsa0FNSI+AOwoCemxUyX0aOxGSNsiQmb1xZvXg455zLmbd5OOecy5kXD+eccznz4uEaHEm9JE2R9JOkeZL+I6lOrqSTtKWk/pLWqca2/SVZ2u0rScOjNrmq9r1fUrY5tp2rE148XIOiMJ3wPcAYoBtwLPAscHAdPeWWQD+gyuIR+RHYLbqdD/wRGC+peRX7XQ70qmFG53LmPcxdQ3M6cJeZXZK2boSkAXEFyrDMzN6Ifn9D0izgFaAz8ETmxpKamtliM5tWyJDO+ZGHa2jWAb7OXGlplx1KKo1OGx0t6aHo9Na3kvpl7iepo6Q3JS2R9I2k2xXmk0ZSe2BEtOnn0WPOyDHvlOhnafSYMyRdL+lSSbOBBdH6VU5bSSqR9KikuZLKJL0v6ei0+5tIGiTpC0k/RwPqdc4xn2ug/MjDNTTvAGdE3+ifr2LMn2uB54HDCGNe9ZM018xuA5C0LTAaGAt0BzYBrgY2JwwH8Q7RuEOEU2RzgJ9zzFsa/UwveEcTRkX4GxX8H5a0AaFXcVmU4QtguyhjypPAzoTTatOAI4DnJO1kZu/mmNM1MF48XENzGvAMcD9gkj4ChhMGlluQse2HZnZy9PuY6AP5Ekl3mNkK4FJgJnCwmS0HkPQD8Jik3czs35I+ifb/TzRcRJUkpf5fbg7cDvwEjMvYrGsVo6WeA7QE2pnZnGjd+LTn2IcwxUB7M5sUrX5R0pZAX+Dw6mR1DZeftnINipm9D2xNaCC/nTA+1aXA5NTppjRPZyw/BWwEbBwt7ww8nSockeHAMmCPGkZcD1ga3T4hFJAeaQUAYHw1htnuCIzO2C/dvoSjmdckNU7dCAVmpxpmdw2IH3m4BsfMfia0RYwAkHQC4QqsE4Cb0zb9NmPX1HJbYFb085uMx14u6Xtg3RrG+5HwwW6ED/evbNVhIL5ZZa9VrQe8Xcn9rYENCUUq0/Is65wrx4uHa/DMbIikQcBWGXdtUMHynLSf5baJxhBaD8g6xlA1LDOzqvprVGdMoe8Jxa0iPwBfAn+pbjDn0vlpK9egRO0WmevWJ7QPZH6jPzRjOdXoPTtafhM4NCoY6ds0Bl6Nln+JfhZ68qjxwAGS2lRy/4bAQjObnHkrXEyXVH7k4RqaD6IhrF8knIYqIVyNVAY8kLHtttFw/cMJV1udAJwVNZYDXEGY9vMZSXcQ2kKuAcaY2b+jbVIN5idLGgaUmdkHdfOnlXMjoQPkK5IGEq622hpobmaDCFeIjQHGSrqGcPXW2oROiU3M7OICZHQJ5sXDNTSXAYcAtxDaJb4GXic0Sn+ese2FQFdC8VhC6MX9r9SdZvahpE7AlYTG9AWEeaQvTNtmpqTzgTOBMwhHLaV18YelM7PvJO0ODAJuAtYEPgWuiu43Sd2AS4CzgU0Jp7LeJcyH7VylfEh25zJE03d+Tpja8/l40zhXnLzNwznnXM68eDjnnMuZn7ZyzjmXMz/ycM45lzMvHs4553LmxcM551zOvHg455zLmRcP55xzOft/KYGeakWS0PUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot exact payoff function (evaluated on the grid of the uncertainty model)\n",
    "x = uncertainty_model.values\n",
    "y = np.maximum(0, x - strike_price)\n",
    "plt.plot(x, y, 'ro-')\n",
    "plt.grid()\n",
    "plt.title('Payoff Function', size=15)\n",
    "plt.xlabel('Spot Price', size=15)\n",
    "plt.ylabel('Payoff', size=15)\n",
    "plt.xticks(x, size=15, rotation=90)\n",
    "plt.yticks(size=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:11.087435Z",
     "start_time": "2020-07-13T23:35:11.082555Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact expected value:\t0.1623\n",
      "exact delta value:   \t0.8098\n"
     ]
    }
   ],
   "source": [
    "# evaluate exact expected value (normalized to the [0, 1] interval)\n",
    "exact_value = np.dot(uncertainty_model.probabilities, y)\n",
    "exact_delta = sum(uncertainty_model.probabilities[x >= strike_price])\n",
    "print('exact expected value:\\t%.4f' % exact_value)\n",
    "print('exact delta value:   \\t%.4f' % exact_delta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluate Expected Payoff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:11.253873Z",
     "start_time": "2020-07-13T23:35:11.247765Z"
    }
   },
   "outputs": [],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "# construct amplitude estimation \n",
    "ae = IterativeAmplitudeEstimation(epsilon=epsilon, alpha=alpha, \n",
    "                                  state_preparation=european_call,\n",
    "                                  objective_qubits=[3], \n",
    "                                  post_processing=european_call_objective.post_processing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:35.384882Z",
     "start_time": "2020-07-13T23:35:11.373858Z"
    }
   },
   "outputs": [],
   "source": [
    "result = ae.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:35.392375Z",
     "start_time": "2020-07-13T23:35:35.387058Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact value:        \t0.1623\n",
      "Estimated value:    \t0.1694\n",
      "Confidence interval:\t[0.1633, 0.1755]\n"
     ]
    }
   ],
   "source": [
    "conf_int = np.array(result['confidence_interval'])\n",
    "print('Exact value:        \\t%.4f' % exact_value)\n",
    "print('Estimated value:    \\t%.4f' % (result['estimation']))\n",
    "print('Confidence interval:\\t[%.4f, %.4f]' % tuple(conf_int))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of constructing these circuits manually, Qiskit's finance module offers the `EuropeanCallExpectedValue` circuit, which already implements this functionality as building block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.finance.applications import EuropeanCallExpectedValue\n",
    "\n",
    "european_call_objective = EuropeanCallExpectedValue(num_uncertainty_qubits,\n",
    "                                                    strike_price,\n",
    "                                                    rescaling_factor=c_approx,\n",
    "                                                    bounds=(low, high))\n",
    "\n",
    "# append the uncertainty model to the front\n",
    "european_call = european_call_objective.compose(uncertainty_model, front=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact value:        \t0.1623\n",
      "Estimated value:    \t0.1708\n",
      "Confidence interval:\t[0.1656, 0.1759]\n"
     ]
    }
   ],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "# construct amplitude estimation \n",
    "ae = IterativeAmplitudeEstimation(epsilon=epsilon, alpha=alpha, \n",
    "                                  state_preparation=european_call,\n",
    "                                  objective_qubits=[3], \n",
    "                                  post_processing=european_call_objective.post_processing)\n",
    "result = ae.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)\n",
    "\n",
    "conf_int = np.array(result['confidence_interval'])\n",
    "print('Exact value:        \\t%.4f' % exact_value)\n",
    "print('Estimated value:    \\t%.4f' % (result['estimation']))\n",
    "print('Confidence interval:\\t[%.4f, %.4f]' % tuple(conf_int))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluate Delta\n",
    "\n",
    "The Delta is a bit simpler to evaluate than the expected payoff.\n",
    "Similarly to the expected payoff, we use a comparator circuit and an ancilla qubit to identify the cases where $S_T > K$.\n",
    "However, since we are only interested in the probability of this condition being true, we can directly use this ancilla qubit as the objective qubit in amplitude estimation without any further approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:35.826269Z",
     "start_time": "2020-07-13T23:35:35.819320Z"
    }
   },
   "outputs": [],
   "source": [
    "from qiskit.finance.applications import EuropeanCallDelta\n",
    "\n",
    "european_call_delta = EuropeanCallDelta(num_uncertainty_qubits, strike_price, bounds=(low, high))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">                                                 \n",
       "state_0: ───────■─────────────────────────────■──\n",
       "                │                             │  \n",
       "state_1: ───────┼────■───────────────────■────┼──\n",
       "         ┌───┐  │    │            ┌───┐  │    │  \n",
       "state_2: ┤ X ├──┼────┼─────────■──┤ X ├──┼────┼──\n",
       "         ├───┤  │    │       ┌─┴─┐└───┘  │    │  \n",
       "state_3: ┤ X ├──┼────┼───────┤ X ├───────┼────┼──\n",
       "         └───┘┌─┴─┐  │       └─┬─┘       │  ┌─┴─┐\n",
       " work_0: ─────┤ X ├──■─────────┼─────────■──┤ X ├\n",
       "              └───┘┌─┴─┐┌───┐  │  ┌───┐┌─┴─┐└───┘\n",
       " work_1: ──────────┤ X ├┤ X ├──■──┤ X ├┤ X ├─────\n",
       "                   └───┘└───┘     └───┘└───┘     </pre>"
      ],
      "text/plain": [
       "                                                 \n",
       "state_0: ───────■─────────────────────────────■──\n",
       "                │                             │  \n",
       "state_1: ───────┼────■───────────────────■────┼──\n",
       "         ┌───┐  │    │            ┌───┐  │    │  \n",
       "state_2: ┤ X ├──┼────┼─────────■──┤ X ├──┼────┼──\n",
       "         ├───┤  │    │       ┌─┴─┐└───┘  │    │  \n",
       "state_3: ┤ X ├──┼────┼───────┤ X ├───────┼────┼──\n",
       "         └───┘┌─┴─┐  │       └─┬─┘       │  ┌─┴─┐\n",
       " work_0: ─────┤ X ├──■─────────┼─────────■──┤ X ├\n",
       "              └───┘┌─┴─┐┌───┐  │  ┌───┐┌─┴─┐└───┘\n",
       " work_1: ──────────┤ X ├┤ X ├──■──┤ X ├┤ X ├─────\n",
       "                   └───┘└───┘     └───┘└───┘     "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "european_call_delta.decompose().draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">     ┌───────┐┌──────┐\n",
       "q_0: ┤0      ├┤0     ├\n",
       "     │       ││      │\n",
       "q_1: ┤1 P(X) ├┤1     ├\n",
       "     │       ││      │\n",
       "q_2: ┤2      ├┤2     ├\n",
       "     └───────┘│  ECD │\n",
       "q_3: ─────────┤3     ├\n",
       "              │      │\n",
       "q_4: ─────────┤4     ├\n",
       "              │      │\n",
       "q_5: ─────────┤5     ├\n",
       "              └──────┘</pre>"
      ],
      "text/plain": [
       "     ┌───────┐┌──────┐\n",
       "q_0: ┤0      ├┤0     ├\n",
       "     │       ││      │\n",
       "q_1: ┤1 P(X) ├┤1     ├\n",
       "     │       ││      │\n",
       "q_2: ┤2      ├┤2     ├\n",
       "     └───────┘│  ECD │\n",
       "q_3: ─────────┤3     ├\n",
       "              │      │\n",
       "q_4: ─────────┤4     ├\n",
       "              │      │\n",
       "q_5: ─────────┤5     ├\n",
       "              └──────┘"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "state_preparation = QuantumCircuit(european_call_delta.num_qubits)\n",
    "state_preparation.append(uncertainty_model, range(uncertainty_model.num_qubits))\n",
    "state_preparation.append(european_call_delta, range(european_call_delta.num_qubits))\n",
    "state_preparation.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:35.843101Z",
     "start_time": "2020-07-13T23:35:35.828358Z"
    }
   },
   "outputs": [],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "# construct amplitude estimation \n",
    "ae_delta = IterativeAmplitudeEstimation(epsilon=epsilon, alpha=alpha, \n",
    "                                        state_preparation=state_preparation,\n",
    "                                        objective_qubits=[num_uncertainty_qubits])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:58.524643Z",
     "start_time": "2020-07-13T23:35:35.848588Z"
    }
   },
   "outputs": [],
   "source": [
    "result_delta = ae_delta.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:58.531941Z",
     "start_time": "2020-07-13T23:35:58.526967Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact delta:    \t0.8098\n",
      "Esimated value: \t0.8068\n",
      "Confidence interval: \t[0.8013, 0.8123]\n"
     ]
    }
   ],
   "source": [
    "conf_int = np.array(result_delta['confidence_interval'])\n",
    "print('Exact delta:    \\t%.4f' % exact_delta)\n",
    "print('Esimated value: \\t%.4f' % result_delta['estimation'])\n",
    "print('Confidence interval: \\t[%.4f, %.4f]' % tuple(conf_int))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:58.967675Z",
     "start_time": "2020-07-13T23:35:58.792732Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td>Qiskit</td><td>None</td></tr><tr><td>Terra</td><td>0.17.0.dev0+4ada179</td></tr><tr><td>Aer</td><td>0.6.1</td></tr><tr><td>Ignis</td><td>0.5.0.dev0+470d8cc</td></tr><tr><td>Aqua</td><td>0.9.0.dev0+5a88b59</td></tr><tr><td>IBM Q Provider</td><td>0.8.0</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.7.7 (default, May  6 2020, 04:59:01) \n",
       "[Clang 4.0.1 (tags/RELEASE_401/final)]</td></tr><tr><td>OS</td><td>Darwin</td></tr><tr><td>CPUs</td><td>2</td></tr><tr><td>Memory (Gb)</td><td>16.0</td></tr><tr><td colspan='2'>Tue Oct 20 10:57:12 2020 CEST</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='width: 100%; background-color:#d5d9e0;padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'><h3>This code is a part of Qiskit</h3><p>&copy; Copyright IBM 2017, 2020.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "%qiskit_version_table\n",
    "%qiskit_copyright"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
